/tomo/pyhst

/*
 * The PyHST program is Copyright (C) 2002-2011 of the
 * European Synchrotron Radiation Facility (ESRF) and
 * Karlsruhe Institute of Technology (KIT).
 *
 * PyHST 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 3 of the License, or
 * (at your option) any later version.
 * 
 * hst 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 this program.  If not, see <http://www.gnu.org/licenses/>.
 */

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

#include <sys/mman.h>


#undef HW_USE_SIMD

#ifdef HW_USE_SIMD
# include <xmmintrin.h>
#endif /* HW_USE_SIMD */

#include "Vhst_calculate_limits.h"

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

#ifdef HW_USE_CUDA
# include "hst_cuda/hst_cuda.h"
# ifdef PYHST_ENABLE_DFI
#  include "dfi_cuda/dfi_cuda.h"
# endif /* PYHST_ENABLE_DFI */
#endif /* HW_USE_CUDA */

#ifdef HW_USE_OPENCL
# include "hst_opencl/hst_opencl.h"
#endif /* HW_USE_OPENCL */

#ifdef HW_USE_CPU
# include "hst_cpu/hst_cpu.h"
#endif /* HW_USE_CPU */

#include "hw_sched.h"

#define PRE_RADIUS 1

static int hst_copies = 0;

#ifdef HW_USE_CUDA
static HSTReconstructor *cuda_recon = NULL;
static HSTReconstructor *cuda_fbp_recon = NULL;
# ifdef PYHST_ENABLE_DFI
static HSTReconstructor *cuda_dfi_recon = NULL;
# endif /* PYHST_ENABLE_DFI */
#endif /* HW_USE_CUDA */

#ifdef HW_USE_OPENCL
static HSTReconstructor *opencl_recon = NULL;
#endif /* HW_USE_OPENCL */

#ifdef HW_USE_CPU
static HSTReconstructor *cpu_recon = NULL;
#endif /* HW_USE_CPU */

int hst_init(void) {
    hw_sched_init();

#ifdef HW_USE_CUDA
    cuda_fbp_recon = hst_cuda_init(HW_USE_CUDA);
# ifdef PYHST_ENABLE_DFI
    cuda_dfi_recon = dfi_cuda_init(HW_USE_CUDA);
# endif /* PYHST_ENABLE_DFI */
#endif /* HW_USE_CUDA */

#ifdef HW_USE_OPENCL
    opencl_recon = hst_opencl_init(HW_USE_OPENCL);
#endif /* HW_USE_OPENCL */

#ifdef HW_USE_CPU
    cpu_recon = hst_cpu_init(HST_RECONSTRUCTOR_USE_CPU);
#endif /* HW_USE_CPU */

    return 0;
}

void hst_free() {
#ifdef HW_USE_CPU
    hst_cpu_free();
#endif /* HW_USE_CPU */

#ifdef HW_USE_CUDA
# ifdef PYHST_ENABLE_DFI
    dfi_cuda_free();
# endif /* PYHST_ENABLE_DFI */
    hst_cuda_free();
#endif /* HW_USE_CUDA */

#ifdef HW_USE_OPENCL
    hst_opencl_free();
#endif /* HW_USE_OPENCL */
}

size_t hst_get_balanced_number_of_reconstructors() {
    size_t count = 0;

#ifdef HW_USE_THREADS
# ifdef HW_USE_CUDA
    count += cuda_fbp_recon->devices;
# endif /* HW_USE_OPENCL */

# ifdef HW_USE_OPENCL
    count += opencl_recon->devices;
# endif /* HW_USE_CUDA */

# if HW_USE_CPU	== 1
    count += cpu_recon->devices - 1; // one core reserved for data preloading
# elif defined(HW_USE_CPU)
    if (!count) count += cpu_recon->devices - 1;
# endif /* HW_USE_CPU */
#else
    count = 1;
#endif

	// we should try to estimate relative performance here and compute LCM
	// this should allow PyHST to allocate such number of slices, that all 
	// reconstructors terminate simultaneously
    return count;
}

/**
 * HST Context
 */
struct HSTContextT {
    HSTSetup setup;			//!< HST paramaters
    
    int num_recon;			//!< Number of initialized reconstructors
    HWSched sched;			//!< Scheduller
    HWSched cpu_sched;                  //!< Preprocessing Scheduller
    HSTReconstructorContext **recon;	//!< Pool of reconstructors (NULL-terminated)
	
	// some cached values set in hst_set_projections
    const float *custom_angles;		//!< angles of projection axis (per projection)
    const float *axis_corrections; 	//!< rotation axis corrections (per projection)

    float *result_tmpbuf;		//!< Temporary buffer for reconstructed slices
    float *result_reqbuf;		//!< Result buffer supplied in reconstruction request
    const float *sinograms;		//!< Sinograms for processing
    HWMutex output_mutex;		//!< Output serializer
    FILE *output;            		//!< Output file (open and close is done in PyHST_c)
#ifdef HST_USE_FASTWRITER
    fastwriter_t *fw;			//!< Output with fastwriter
#else /* HST_USE_FASTWRITER */
    void *fw;
#endif /* HST_USE_FASTWRITER */

    int limits_configured;		//!< Flag indicating if limits are already configured
    int filter_configured;		//!< Flag indicating if filter is already configured

    int num_projections;
    int current_projection;
    int preproc_slices;
    float *preproc_sinograms;
    float **projection_data;

    int num_slices;
    int current_slice;

    int slice_offset;			//!< Number of slices already processed (for navigation in output file)
    int processing;			//!< Flag indicating the working threads are running and processing slices

    void *mmap;

#ifdef PYHST_MEASURE_TIMINGS
    double comp_timer;
    double io_timer;
#endif /* PYHST_MEASURE_TIMINGS */

};
typedef struct HSTContextT HSTContext; //!< a short name


static int hst_configure_data_structures(HSTContext *ctx, const float *custom_angles, const float *axis_corrections) {
    int projection;		// projection counter
    int num_projections;	// Number of projections

    float angle;		// projections angle
    float angle_offset;         // Offset rotation angle in radians, could be overriden with custom_angles
    float angle_increment;      // Angle in radians between each, could be overriden with custom_angles
    float *cos_s, *sin_s;

    float axis_position;       	// Position of rotation axis
    float *axis_position_corr_s;

    ctx->custom_angles = custom_angles;
    ctx->axis_corrections = axis_corrections;

    num_projections = ctx->setup.num_projections;

    angle_offset = ctx->setup.angle_offset;
    angle_increment = ctx->setup.angle_increment;
    axis_position = ctx->setup.axis_position;

	// Computing FFT dimensions, hard to understand why the -1 is there but it works nicely
	// DS: This is actually two times more than the closest power of 2. Do we need that?
    ctx->setup.dim_fft = 1<<((int)(log((1. * ctx->setup.fft_oversampling * ctx->setup.num_bins - 1)) / log(2.0) + 0.9999));
    printf("FFT convolution: %u\n", ctx->setup.dim_fft);
    ctx->setup.filter = malloc((2 * ctx->setup.dim_fft) * sizeof(float));
    check_alloc(ctx->setup.filter, " FILTER");

    axis_position_corr_s = (float *) malloc(num_projections * sizeof(float));
    check_alloc(axis_position_corr_s, "axis_position_corr_s");
    ctx->setup.axis_position_corr_s = axis_position_corr_s;
    
    cos_s = (float *) malloc(num_projections * sizeof(float));
    check_alloc(cos_s, "cos_s");
    ctx->setup.cos_s = cos_s;

    sin_s = (float *) malloc(num_projections * sizeof(float));
    check_alloc(sin_s, "sin_s");
    ctx->setup.sin_s = sin_s;

    for (projection = 0; projection < num_projections; projection++) {
        if (custom_angles) {
            angle = custom_angles[projection];
        } else {
            angle = angle_offset + projection * angle_increment ;
        }

        cos_s[projection] = cos(angle);
        sin_s[projection] = sin(angle);

	if (axis_corrections) axis_position_corr_s[projection] = axis_position + axis_corrections[projection];
	else axis_position_corr_s[projection] = axis_position;
    }

    return 0;
}

static int hst_configure_limits(HSTContext *ctx) {
    int i;
    
    int projection;		// projection counter
    
	// Some variables from configuration
    int num_projections;	// Number of projections
    int num_bins;		// Number of bins in sinograms

    int start_x;                // First X-pixel in region required by user
    int start_y;                // First Y-pixel in region required by user
    int num_x;			// Number of X-pixels in reconstruction region
    int num_y;			// Number of Y-pixels in reconstruction region

    float axis_position;       	// Position of rotation axis
    float angle_increment;      // Angle in radians between each, could be overriden with custom_angles

    int axis_to_the_center;
    
    int zerooffmask;

    int *minX = NULL, *maxX = NULL;
    float *cos_s, *sin_s;
    float *axis_position_corr_s;

	// The block of variables for reconstruction limits
    int status;
    int *X_STARTS, *X_ENDS;
    int y_start;
    int y_end;
    int startxdum, startydum, finxdum, finydum;

	// The next block are variables needed for zeroofmask-mode computations 
#if PRE_RADIUS
    float * RRadius = NULL;
#endif /* PRE_RADIUS */
    float yy, rradius;
    int yoristart, yoriend;
    int yorigin ;
    int x_start, x_end ;
    int y;
    float EXT1, EXT2;
    int XH;


    /* The block of variables describing the position of axes (and used to move
    the axis to the center of the reconstructed zone) as well as projection
    angles */
    float MOVE_X;
    float MOVE_Y;

    int fai360;
    float axis_position_corr = 0;
    float cos_angle , sin_angle;
    float corre; // axis correction

    assert(ctx);
    
    num_bins = ctx->setup.num_bins;
    num_projections = ctx->setup.num_projections;
    
    start_x = ctx->setup.start_x;
    num_x = ctx->setup.num_x;
    start_y = ctx->setup.start_y;
    num_y = ctx->setup.num_y;

    axis_position_corr_s = ctx->setup.axis_position_corr_s;
    cos_s = ctx->setup.cos_s;
    sin_s = ctx->setup.sin_s;
    
    axis_position = ctx->setup.axis_position;
    angle_increment = ctx->setup.angle_increment;

    axis_to_the_center = ctx->setup.center_axis;
    zerooffmask = ctx->setup.zerooffmask;

    
    /*************************************************************
     * Depending on geometry ( dimensions, axis positions ....  )*
     * only a restraint part of the slice can be back projected  *
     * from all the sinograms.                                   *
     * X_Start tells the limits for each y  ( idem for X_ENDS)   *
     *************************************************************/
    X_STARTS = malloc(num_bins * sizeof(int));
    check_alloc(X_STARTS," X_STARTS allocation in c_hst_recon_overslices");
    X_ENDS = malloc(num_bins * sizeof(int));
    check_alloc(X_ENDS," X_ENDS allocation in c_hst_recon_overslices");

    /************************************************************
     * calculate limits of slices to reconstruct                *
     ************************************************************/
    /* Restrict the area if projections are covering more than 360 grads  */
    startxdum=1;
    startydum=1;
    finxdum = num_bins ;
    finydum = num_bins ;

    /* DSBug: case of custom angles is not handled here */
    if (fabs((num_projections)*(angle_increment) - 2*M_PI)>0.001 ) {
        hst_calculate_limits__(&num_bins, &startxdum, &startydum, &finxdum, &finydum, &axis_position, &y_start, &y_end, X_STARTS, X_ENDS, &status);
        y_start+=1;
        y_end  -=1;
    } else {
        y_start=-1;
        y_end  =-1;
    }

    fai360 = (y_start == -1 && y_end == -1);
    ctx->setup.fai360 = fai360;


    if (axis_to_the_center) {
        MOVE_X = (((start_x - 1) + num_x / 2.0) - axis_position);
        MOVE_Y = (((start_y - 1) + num_y / 2.0) - axis_position);
    } else {
        MOVE_X = 0;
        MOVE_Y = 0;
    }

    ctx->setup.offset_x = start_x - MOVE_X ;
    ctx->setup.offset_y = start_y - MOVE_Y ;

    if (zerooffmask) {
        /* these should be allocated once and for all out of the routine */
	if (ctx->setup.minX) minX = ctx->setup.minX;
	else {
    	    minX = (int*) malloc(num_y * sizeof(int));
    	    check_alloc(minX, "minX");
	    ctx->setup.minX = minX;
	}
	
	if (ctx->setup.maxX) maxX = ctx->setup.maxX;
	else {
    	    maxX = (int*) malloc(num_y * sizeof(int));
    	    check_alloc(maxX, "maxX");
	    ctx->setup.maxX = maxX;
	}
	
        for (i = 0; i < num_y; i++) {
            maxX[i] = 0;
            minX[i] = 1000000;
        }

#if PRE_RADIUS
        /* The difference with !pre_radius what we do not suing axis_corrections,
        i.e axis_position is used in place of axis_position_corr */
        RRadius = (float *) malloc(num_y * sizeof(float));
        check_alloc(RRadius, "RRadius");

        for (y = 0; y < num_y; y++) {
            rradius = MAX(axis_position , num_bins - axis_position) ;
            yy = ((y + start_y - MOVE_Y - 1) + 0.5 - axis_position);
            if (fabs(yy) > rradius) {
                RRadius[y] = 0.0;
            } else {
                RRadius[y] = sqrt(rradius * rradius - yy * yy);
            }
        }
#endif /* PRE_RADIUS */
    }

    for (projection = 0; projection < num_projections; projection++) {
        cos_angle = cos_s[projection];
        sin_angle = sin_s[projection];
	axis_position_corr = axis_position_corr_s[projection];
	
	//corre = axis_position_corr - axis_position;
	corre = ctx->axis_corrections?ctx->axis_corrections[projection]:0;

        if (zerooffmask) {
            if (fai360) {
                yoristart = 0;
                yoriend = num_y ;
            } else {
                yoristart = y_start - 1 ;
                yoriend = y_end ;
            }

            for (yorigin = yoristart ; yorigin < yoriend ; yorigin++) {
                if (fai360) {
                    y = yorigin;
                } else {
                    /* y_start e y_end sono calcolati per una maschera massimale */
                    /* y e relativo alla finestra dell user */
                    y = (yorigin + MOVE_Y + corre - start_y) + 100000; /* per evitare che -0.1 sia tagliato a 0 come 0.9 */
                    y = y - 100000;
                }

                /* si potrebbe eventualmente togliere il 2 se tutto va bene : ps 2 ->1*/
                if (y < 0 || y >= num_y - 0) continue;

                if (fai360) {
#if PRE_RADIUS
                    rradius = RRadius[y];
#else /*PRE_RADIUS */
                    rradius = MAX(axis_position_corr , num_bins - axis_position_corr) ;
                    yy = ((y + start_y - MOVE_Y - 1) + 0.5 - axis_position);
                    if (fabs(yy) > rradius) {
                        rradius = 0.0;
                    } else {
                        rradius = sqrt(rradius * rradius - yy * yy);
                    }
#endif /* PRE_RADIUS */


                    if (fabs(cos_angle) > 1.0e-6) { /* s puo lasciare solo la parte in y */
                        EXT1 = -((float)(((axis_position_corr +
                                           ((start_x - MOVE_X - 1) + 0.5 - axis_position) * cos_angle -
                                           ((y + start_y - MOVE_Y - 1) + 0.5 - axis_position) * sin_angle - 0.5)))) / cos_angle ;
                        EXT2 = (num_bins - (float)(((axis_position_corr +
                                                     ((start_x - MOVE_X - 1) + 0.5 - axis_position) * cos_angle -
                                                     ((y + start_y - MOVE_Y - 1) + 0.5 - axis_position) * sin_angle - 0.5)))) / cos_angle ;
                    } else {
                        XH = ((float)(((axis_position_corr -
                                        ((y + start_y - MOVE_Y - 1) + 0.5 - axis_position) * sin_angle - 0.5)))
                             );
                        if (XH <= 0 || XH >= num_bins) continue;
                        EXT1 = 0;
                        EXT2 = num_x;
                    }
                    x_start = MIN(EXT1, EXT2) + 1;
                    x_end = MAX(EXT1, EXT2) - 1;

                    /********************************************************************************************/
                    x_start = MAX(x_start , axis_position - rradius - start_x); /* SOSPETTO */
                    x_end = MIN(x_end , axis_position + rradius - start_x);
                    /********************************************************************************************/

                } else { // fai360 if
                    x_start = (X_STARTS[yorigin] - 1) - start_x + MOVE_X + corre ;
                    /* if( x_start< + DX0) x_start= DX0; */
                    x_end = (X_ENDS[yorigin] - 1) - start_x + MOVE_X + corre ;
                } // end of fai360 if

                {
                    if (x_start < 0) x_start = 0;
                    if (x_end >= num_x - 0) x_end = num_x - 0;
                    /* anche qui se tutto va bene si puo provare ad essere piu larghi
                    x_start+=2;
                    x_end-=2; */
                }

                if (zerooffmask) {
                    minX[y] = MIN(minX[y], x_start);
                    maxX[y] = MAX(maxX[y], x_end);
                }
            } // end of yorigin loop
        } // zerooffmask if
    } // end of projections loop

#if PRE_RADIUS
    if (zerooffmask) {
	free(RRadius);
    }
#endif /* PRE_RADIUS */

    free(X_ENDS);
    free(X_STARTS);

    ctx->limits_configured = 1;
    ctx->setup.what_changed |= HST_LIMITS_CHANGED;
    
    return 0;
}

#include "astra_fourier.h"
int hst_set_filter(HSTContext *ctx, HSTFilterFunction func, void *func_param) {
    int j;
    int dim_fft;
    float *FILTER;

    assert(ctx);

    dim_fft = ctx->setup.dim_fft;
    FILTER = ctx->setup.filter;

#ifdef PYHST_ASTRA_SCALING
	// dim_fft is next power of two of (2 * num_proj)
	// num_bins is dim_fft/2 + 1
    int num_bins = dim_fft / 2 + 1;	// number of non-symmetric coefficients

    memset(FILTER, 0, 2 * dim_fft * sizeof(float));
    FILTER[0] = 0.25;
    FILTER[1] = 0;
    for (j = 1; j < dim_fft; j++) {
	if (j&1) {
	    double coef;
	    if (2 * j > dim_fft)
		coef = M_PI * (dim_fft - j);
	    else
		coef = M_PI * j;
	    FILTER[2 * j] = -1. / (coef * coef);
	} else {
	    FILTER[2 * j] = 0;
	}
	FILTER[2 * j + 1] = 0;
    }

    int *ip = (int*)malloc((2 + sqrt(dim_fft) + 1)*sizeof(int)); ip[0] = 0;
    float *w = (float*)malloc(dim_fft * sizeof(float) / 2);
    cdft(2 * dim_fft, -1, FILTER, ip, w);
    free(w);
    free(ip);

    for (j = 0; j < num_bins; j++) {
	float _fD = 1.0f;
        float fRelIndex = (float)j / (float)dim_fft;
        float pfW = 2. * M_PI * fRelIndex;

	    // On last iteration it is actually equal. So, rounding error works differently with ASTRA and here.... Therefore, we require a bit on top.
	if (pfW > (M_PI * _fD + 1E-6)) {
	    FILTER[2 * j] = 0;
	    FILTER[2 * j + 1] = 0;
	} else {
	    FILTER[2 * j] = 2.0 * FILTER[2 * j];
	    FILTER[2 * j + 1] = FILTER[2 * j];
	}
    }

    for (j = num_bins; j < dim_fft; j++) {
	FILTER[2 * j] = 0;
	FILTER[2 * j + 1] = 0;
    }

/*
    for (j = 0; j < num_bins; j++) {
	printf("%8.4f %8.4f ", FILTER[2 * j], FILTER[2 * j + 1]);
    }
    printf("\n");
*/
#else /* ASTRA_SCALING */
    if (func) {
        func(func_param, dim_fft, FILTER);
    } else {
        FILTER[0]=0.0;
        FILTER[1]= 1.0/dim_fft ;             // This is the real part of the highest frequency 

        for (j=1;j<dim_fft / 2; ++j) {
            FILTER[2 * j] = j * 2.0 / dim_fft / dim_fft;
            FILTER[2 * j + 1] = FILTER[2 * j];
        }
    }

    if (ctx->setup.sum_rule == 2) {
        FILTER[0]=FILTER[2]/4.0;
    }

    if (ctx->setup.no_filtering) {
        FILTER[1] = -9999.9999;
    }
    /* This is due to different organization of data representation used after by
    real to complex FFT transfer. In both cases only non-redundant coefficients
    are stored however, the current CPU code stores High Frequency Term (HFT) in the
    res[1], while cuFFT (and FFTW as well) have always res[1] = 0 and stores
    the HFT in res[dim_fft] */

    FILTER[dim_fft] = FILTER[1]; // or should it be 0?

    /* CUDA code now includes optimization allowing to compute two convolutions 
    parallely using single Complex-Complex Fourier transform. For that we need
    to feel complete FILTER matrix including redundant coefficients. Therfore, 
    the filter is filled withcoefficients according to the cuFFT/FFTW approach. 
    The case of embedded FFT code (when compressed FFT represntation is needed)
    is handled separately in hst.c/hst_cpu code. */

    FILTER[dim_fft + 1] = FILTER[1];
    FILTER[1] = FILTER[0];
#endif /* ASTRA_SCALING */

    for (j = dim_fft + 2; j < 2 * dim_fft; j+=2) {
	FILTER[j] = FILTER[2 * dim_fft - j];
	FILTER[j + 1] = FILTER[2 * dim_fft - j + 1];
    }


    ctx->filter_configured = 1;
    ctx->setup.what_changed |= HST_FILTER_CHANGED;
    
    return 0;
}

int hst_set_output_file(HSTContext *ctx, FILE *output) {
    if ((ctx->fw)||(output != ctx->output)) {
	ctx->fw = NULL;
	ctx->output = output;
	ctx->setup.what_changed |= HST_OUTPUT_CHANGED;
    }
    
    return 0;
}

#ifdef HST_USE_FASTWRITER
int hst_set_fastwriter(HSTContext *ctx, fastwriter_t *output) {
# ifndef HST_MULTISLICE_MODE
    pyhst_error("Fastwriter is not supported in non-multislice mode");
    return -1;
# endif /* HST_MULTISLICE_MODE */

    if ((ctx->output)||(ctx->fw != output)) {
	ctx->output = NULL;
	ctx->fw = output;
	ctx->setup.what_changed |= HST_OUTPUT_CHANGED;
    }
}
#endif /* HST_USE_FASTWRITER */


int hst_set_padding(HSTContext *ctx, int zero_padding) {
    if (zero_padding != ctx->setup.zero_padding) {
	ctx->setup.zero_padding = zero_padding;
	ctx->setup.what_changed |= HST_PADDING_CHANGED;
    }
    return 0;
}

int hst_set_axis_mode(HSTContext *ctx, int axis_to_the_center) {
    if (axis_to_the_center != ctx->setup.center_axis) {
        ctx->setup.center_axis = axis_to_the_center;
	    // The projections configuration is dependent on this parameter
	hst_configure_limits(ctx);
    }
    return 0;
}

static int hst_threads_configure(HWThread thr, void *hwctx, int device_id, void *data) {
    int err = 0;
    HSTContext *ctx = (HSTContext*)data;

    //printf("Configure %i\n", device_id);

    HSTReconstructorContextPtr rctx = ctx->recon[device_id];
    
    if (rctx->recon.configure) {
	err = rctx->recon.configure(rctx, ctx->setup.what_changed);
    }

    return err;
}

static int hst_configure(HSTContext *ctx) {
    int err;

    assert(ctx);

    if (!ctx->filter_configured) {
        err = hst_set_filter(ctx, NULL, NULL);
	check_code(err, "hst_set_filter");
    }

    if (!ctx->limits_configured) {    
        err = hst_configure_limits(ctx);
	check_code(err, "hst_configure_limits");
    }

    err = hw_sched_execute_thread_task(ctx->sched, ctx, hst_threads_configure);
    check_code(err, "hst_configure");
    
/*
    int i;
    for (i = 0; i < ctx->num_recon; ++i) {
	rctx = ctx->recon[i];
	if (rctx->recon.configure) {
	    err = rctx->recon.configure(rctx, ctx->setup.what_changed);
	    check_code(err, "hst_configure");
	}
    }
*/

    if (ctx->setup.what_changed&HST_OUTPUT_CHANGED) {
	ctx->slice_offset = 0;
    }
    
    ctx->setup.what_changed = 0;
    
    return 0;
}

static int hst_reconstruct_slice(HSTContext *ctx, HSTReconstructorContextPtr rctx, float *slice_result, const float *sinograms ) {
    int err;

    assert(ctx);
    assert(rctx);
    assert(rctx->recon.reconstruct);

#ifdef HST_MULTISLICE_MODE
        if (rctx->recon.send) {
            const float *last_sino;
            float *last_slice;

            last_slice = rctx->saved_slice;
            last_sino = rctx->saved_sino;

            if (sinograms) 
                rctx->recon.send(rctx, sinograms);

            rctx->saved_slice = slice_result;
            rctx->saved_sino = sinograms;

            slice_result = last_slice;
            sinograms = last_sino;
        } 

        if (!slice_result) return 0;
#endif /* HST_MULTISLICE_MODE */
//

    if (rctx->recon.preprocess) {
	err = rctx->recon.preprocess(rctx, slice_result, sinograms);
	if (err) return err;
    }

    err = rctx->recon.reconstruct(rctx, slice_result, sinograms);
    if (err) return err;

    if (rctx->recon.postprocess) {
	err = rctx->recon.postprocess(rctx, slice_result, sinograms);
	if (err) return err;
    }

    rctx->processed_slices += SLICE_BLOCK;

    return 0;
}

static int hst_threads_reconstruct_slice(HWThread thr, void *hwctx, int slice_id, void *data) {
    int err = 0;
    int device_id = thr->thread_id;
    HSTContext *ctx = (HSTContext*)data;
    int cnt, count =  ctx->setup.num_x * ctx->setup.num_y;

#ifdef PYHST_MEASURE_TIMINGS
    GTimer *timer;
#endif /* PYHST_MEASURE_TIMINGS */

#ifndef PYHST_IO_BENCHMARK
    if (slice_id == HW_SCHED_CHUNK_FREE) {
        err = hst_reconstruct_slice(ctx, ctx->recon[device_id], NULL, NULL);
        if (ctx->recon[device_id]->recon.wait)
            ctx->recon[device_id]->recon.wait(ctx->recon[device_id]);
    } else if (slice_id >= 0) {
        err = hst_reconstruct_slice(ctx, ctx->recon[device_id],
# ifdef HST_MULTISLICE_MODE
	    (ctx->result_reqbuf?ctx->result_reqbuf:ctx->result_tmpbuf) + SLICE_BLOCK * slice_id * count,
# else /* HST_MULTISLICE_MODE */
	    ctx->result_reqbuf?(ctx->result_reqbuf + SLICE_BLOCK * slice_id * count):(ctx->result_tmpbuf + SLICE_BLOCK * device_id * count),
# endif /* HST_MULTISLICE_MODE */
	    ctx->sinograms + SLICE_BLOCK * slice_id * ctx->setup.num_bins * ctx->setup.num_projections
        );
    }
    if (err) return err;
#endif /* ! PYHST_IO_BENCHMARK */


#ifndef PYHST_RECON_BENCHMARK
    if (!ctx->result_reqbuf) {
	hw_sched_lock_mutex(ctx->output_mutex);

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

# ifdef HST_MULTISLICE_MODE
        if (slice_id == HW_SCHED_CHUNK_TERMINATOR) {
#  ifdef HST_USE_FASTWRITER  
    	    if (ctx->fw) {
retry:
    		err = fastwriter_push_data(ctx->fw, SLICE_BLOCK * ctx->num_slices * count * sizeof(float), ctx->result_tmpbuf);
    		if (err) {
    		    if (err == EWOULDBLOCK) goto retry;
#   ifdef PYHST_MEASURE_TIMINGS
	            if (timer) g_timer_destroy(timer);
#   endif
	            hw_sched_unlock_mutex(ctx->output_mutex);
	            pyhst_error("faswriter failed wrting %u MB of data, error %i", SLICE_BLOCK * ctx->num_slices * count * sizeof(float)/1024/1024, err);
	            return 1;
    		}
    	    } else {
#  endif /* HST_USE_FASTWRITER */
        	cnt = fwrite(ctx->result_tmpbuf, sizeof(float), SLICE_BLOCK * ctx->num_slices * count, ctx->output);
                if (cnt !=  SLICE_BLOCK * ctx->num_slices * count) {
#  ifdef PYHST_MEASURE_TIMINGS
	            if (timer) g_timer_destroy(timer);
#  endif
	            hw_sched_unlock_mutex(ctx->output_mutex);
	            pyhst_error("fwrite failed, send pushed characters %i, but written %i, errno %i", SLICE_BLOCK * ctx->num_slices * count, cnt, errno);
	            return 1;
	        }
#  ifdef HST_USE_FASTWRITER  
            }
#  endif /* HST_USE_FASTWRITER */
        }
# else /* HST_MULTISLICE_MODE */
	err = fseek(ctx->output, ((long)(ctx->slice_offset + slice_id * SLICE_BLOCK)) * count * sizeof(float), SEEK_SET);
	if (err) {
#  ifdef PYHST_MEASURE_TIMINGS
	    if (timer) g_timer_destroy(timer);
#  endif
	    hw_sched_unlock_mutex(ctx->output_mutex);
	    pyhst_error("fseek failed, errno: %i", errno);
	    return err;
	}

        cnt = fwrite(ctx->result_tmpbuf + device_id * count, sizeof(float), SLICE_BLOCK * count, ctx->output);
	if (cnt != SLICE_BLOCK * count) {
#  ifdef PYHST_MEASURE_TIMINGS
	    if (timer) g_timer_destroy(timer);
#  endif
	    hw_sched_unlock_mutex(ctx->output_mutex);
	    pyhst_error("fwrite failed, send pushed characters %i, but written %i, errno %i", SLICE_BLOCK * count, cnt, errno);
	    return 1;
	}
# endif /* HST_MULTISLICE_MODE */
	//fsync(fileno(ctx->output));


# ifdef PYHST_MEASURE_TIMINGS
	if (timer) {
		// We should be inside the lock, otherwise multiple accesses.
	    g_timer_stop(timer);
	    ctx->io_timer += g_timer_elapsed(timer, NULL);
	    //ctx->comp_timer -= g_timer_elapsed(timer, NULL);
	    g_timer_destroy(timer);
	}
# endif /* PYHST_MEASURE_TIMINGS */

	hw_sched_unlock_mutex(ctx->output_mutex);
    }
/*    
    FILE *f;
    char tmp[64];
    sprintf(tmp, "/tmp/test/%i-%i",slice_id + ctx->slice_offset, device_id);
    f=fopen(tmp, "w");
    fwrite(ctx->result_tmpbuf + device_id * count, sizeof(float), count, f);
    fclose(f);
*/    
#endif /* ! PYHST_RECON_BENCHMARK */

    return 0;
}


int hst_reconstruct_slices(HSTContextPtr ctx, int num_slices, float *result, const float *sinograms) {
    int err;
#ifdef PYHST_MEASURE_TIMINGS
    GTimer *timer;
#endif /* PYHST_MEASURE_TIMINGS */

    assert(ctx);

    if (ctx->setup.what_changed) {
	hst_configure(ctx);
    }

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

    ctx->num_slices = num_slices / SLICE_BLOCK;
    ctx->result_reqbuf = result;
    ctx->sinograms = sinograms;

    pyhst_info("num_projection: %d, angle_increment: %e, fai360: %s", ctx->setup.num_projections, ctx->setup.angle_increment, ctx->setup.fai360?"on":"off" );
    hw_sched_schedule_task(ctx->sched, ctx, hst_threads_reconstruct_slice);
    err = hw_sched_wait_task(ctx->sched);
    ctx->slice_offset += SLICE_BLOCK * num_slices;

#ifdef PYHST_MEASURE_TIMINGS
    if (timer) {
	ctx->comp_timer += g_timer_elapsed(timer, NULL);
	g_timer_destroy(timer);
    }
#endif /* PYHST_MEASURE_TIMINGS */

/*
    int i;
    for (i=0; i<num_slices; i++) {
	err = hst_reconstruct_slice(ctx,
	    ctx->recon[0],
	    result?result:ctx->result_tmpbuf,
	    sinograms + i * ctx->setup.num_bins * ctx->setup.num_projections
	);
	if (err) return err;
	
        if (!result) {
	    fwrite(ctx->result_tmpbuf, ctx->setup.num_x * sizeof(float), ctx->setup.num_y, ctx->output);
	}

    }
*/    

    return err;
}

int hst_reconstruct_start(HSTContextPtr ctx, int num_slices, const float *sinograms) {
#ifdef PYHST_MEASURE_TIMINGS
    GTimer *timer;
#endif /* PYHST_MEASURE_TIMINGS */

    assert(ctx);

    if (ctx->setup.what_changed) {
	hst_configure(ctx);
    }

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

    ctx->num_slices = num_slices / SLICE_BLOCK;
    ctx->result_reqbuf = NULL;
    ctx->sinograms = sinograms;
        
    pyhst_info("num_projection: %d, angle_increment: %e, fai360: %s", ctx->setup.num_projections, ctx->setup.angle_increment, ctx->setup.fai360?"on":"off" );
    hw_sched_schedule_task(ctx->sched, ctx, hst_threads_reconstruct_slice);
    ctx->processing = 1;

#ifdef PYHST_MEASURE_TIMINGS
    if (timer) {
	ctx->comp_timer += g_timer_elapsed(timer, NULL);
	g_timer_destroy(timer);
    }
#endif /* PYHST_MEASURE_TIMINGS */

#ifndef HW_USE_PARALLEL_IO
    hst_reconstruct_wait(ctx);
#endif /* ! HW_USE_PARALLEL_IO */
    
    return 0;
}

int hst_reconstruct_wait(HSTContextPtr ctx) {
    int err = 0;

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

    if (ctx->processing) {
	err = hw_sched_wait_task(ctx->sched);
	ctx->slice_offset += SLICE_BLOCK * ctx->num_slices;
	ctx->processing = 0;
    }

#ifdef PYHST_MEASURE_TIMINGS
    if (timer) {
	ctx->comp_timer += g_timer_elapsed(timer, NULL);
	g_timer_destroy(timer);
    }
#endif /* PYHST_MEASURE_TIMINGS */

    return err;
}


#undef HW_USE_SIMD
static int hst_threads_preprocess_projection(HWThread thr, void *hwctx, int proj, void *data) {
    HSTContext *ctx = (HSTContext*)data;

    int i,j;
    int num_slices = ctx->preproc_slices;
    int num_bins = ctx->setup.num_bins;
    int num_proj = ctx->setup.num_projections;
    float norm = ctx->setup.norm;
    
    int num_bins4 = (num_bins/4);

#ifdef HW_USE_SIMD
    __m128 norm4 = _mm_set_ps1(norm);
#endif /* HW_USE_SIMD */
    
    while (!ctx->projection_data[proj]) 
        usleep(100);

//    printf("[%p] %i, %i of %i, %i %f\n", thr, num_slices, proj, num_proj, num_bins, norm);
    for (i = 0; i < num_slices; i++) {
        int src_offset = i * num_bins;
        int dst_offset = i * num_bins * num_proj + proj * num_bins;

#ifdef HW_USE_SIMD
            // Looks pretty inoptimal... Enforce padding? Hand-written assembler?
        __m128* src = (__m128*)(ctx->projection_data[proj] + src_offset);
        __m128* dst = (__m128*)(ctx->preproc_sinograms + dst_offset);

        for (j = 0; j < num_bins4; j++) {
            dst[j] = _mm_mul_ps(src[j], norm4);
        }

        for (j = num_bins4 * 4; j < num_bins; j++) {
#else /* HW_USE_SIMD */
        for (j = 0; j < num_bins; j++) {
#endif /* HW_USE_SIMD */
            ctx->preproc_sinograms[dst_offset + j] = norm * ctx->projection_data[proj][src_offset + j];
        }
    }

    return 0;
}

int hst_preprocess_start(HSTContextPtr ctx, int num_slices, const float *sinograms) {
    int err = 0;
    
    ctx->num_projections = ctx->setup.num_projections;
    ctx->preproc_slices = num_slices;
    ctx->preproc_sinograms = sinograms;

#ifdef HW_USE_PARALLEL_IO
    memset(ctx->projection_data, 0, ctx->setup.num_projections * sizeof(float*));
    err = hw_sched_schedule_task(ctx->cpu_sched, ctx, hst_threads_preprocess_projection);
#endif /* ! HW_USE_PARALLEL_IO */

    return err;
}

int hst_preprocess_wait(HSTContextPtr ctx) {
    int err = 0;
#ifdef HW_USE_PARALLEL_IO
    err = hw_sched_wait_task(ctx->cpu_sched);
#endif /* ! HW_USE_PARALLEL_IO */
    return err;
}

int hst_preprocess_projection(HSTContextPtr ctx, int proj, float *data) {

    ctx->projection_data[proj] = data;

#ifndef HW_USE_PARALLEL_IO
    hst_threads_preprocess_projection(NULL, NULL, proj, ctx);
#endif /* ! HW_USE_PARALLEL_IO */

/*
  int i,j;
  int num_slices = ctx->preproc_slices;
  int num_bins = ctx->setup.num_bins;
  int num_proj = ctx->setup.num_projections;
  float norm = ctx->setup.norm;

//  printf("%i, %i of %i, %i %f\n", num_slices, proj, num_proj, num_bins, norm);
  for (i = 0; i < num_slices; i++) {
    for (j = 0; j < num_bins; j++) {
        ctx->preproc_sinograms[i * num_bins * num_proj + proj * num_bins + j] = norm * data[i * num_bins + j];
    }
  }
*/
  return 0;
}


HSTReconstructorConstContextPtr *hst_get_configured_reconstructors(HSTContextPtr ctx) {
    assert(ctx);
    
    return (HSTReconstructorConstContextPtr*)ctx->recon;
}

#ifdef PYHST_MEASURE_TIMINGS
double hst_get_io_timer(HSTContextPtr ctx) {
    assert(ctx);
    return ctx->io_timer;
}

double hst_get_recon_timer(HSTContextPtr ctx) {
    assert(ctx);
    return ctx->comp_timer;
}
#endif /* PYHST_MEASURE_TIMINGS */


HSTContext *hst_create_context(void) {
    if (hst_copies) pyhst_fail("Only a single HST context is supported by the current version");
    
    HSTContext *ctx = (HSTContext*)malloc(sizeof(HSTContext));
    if (ctx) {
        memset(ctx, 0, sizeof(HSTContext));
	++hst_copies;
    }
    return ctx;
}

static int hst_threads_init(HWThread thr, void *hwctx, int device_id, void *data) {
    int err;
    
    HSTContext *ctx = (HSTContext*)data;

    //printf("Init %i\n", device_id);
    err = hst_reconstructor_init(ctx->recon[device_id], thr);
    check_code(err, "hst_reconstructor_init");
    
    return 0;
}

void *hst_pinned_malloc(size_t size, HSTPinnedAccess a) {
#ifdef HW_USE_OPENCL
    return hst_opencl_host_malloc(size, a);
#else /* OPENCL */
# ifdef HW_USE_CUDA
    return hst_cuda_host_malloc(size, a);
# else /* USEGPU */
    return malloc(size);
/*
    Unfortunatelly, this demands ROOT priveleges (or non-standard system
    configuration) and anyway not supported by CUDA 2.2, so there is only
    way to allocate through CUDA (probably we should get results in our
    personal buffers to avoid such problems.
    if (ctx->result_tmpbuf) mlock(ctx->result, 1 * ctx->setup.num_x * ctx->setup.num_y * sizeof(float));
    ctx->result_tmpbuf = mmap(NULL, 1 * ctx->setup.num_x * ctx->setup.num_y * sizeof(float), PROT_READ|PROT_WRITE, MAP_SHARED|MAP_ANONYMOUS|MAP_LOCKED, -1, 0);
*/
# endif /* USEGPU */
#endif /* OPENCL */
}

void hst_pinned_free(void *ptr) {
#ifdef HW_USE_OPENCL
    hst_opencl_host_free(ptr);
#else /* OPENCL */
# ifdef HW_USE_CUDA
    hst_cuda_host_free(ptr);
# else /* USEGPU */
    free(ptr);
# endif /* USEGPU */
#endif /* OPENCL */
}

int hst_init_context(HSTContext *ctx, HSTSetup *setup, const float *custom_angles, const float *axis_corrections) {
    int err;
    int i, limit, threads;
    size_t size;

    assert(ctx);
    assert(setup);

    hst_setup_init(&ctx->setup);
    
    size = ((void*)&setup->END_OF_GLOBAL_PART) - ((void*)setup);
    memcpy(&ctx->setup, setup,  size);

    ctx->setup.what_changed = HST_ALL_CHANGED;
    hst_configure_data_structures(ctx, custom_angles, axis_corrections);


# ifdef HW_USE_CUDA
    switch (setup->method) {
     case HST_METHOD_FBP: 
        cuda_recon = cuda_fbp_recon;
        break;
# ifdef PYHST_ENABLE_DFI
     case HST_METHOD_DFI:
        cuda_recon = cuda_dfi_recon;
        break;
# endif /* PYHST_ENABLE_DFI */
      default:
        pyhst_error("Unknown reconstruction mode specififed");
        exit(-1);
    }
#endif /* HW_USE_CUDA */

#ifdef HW_USE_THREADS
    threads = 0;

# ifdef HW_USE_CUDA
    threads += cuda_recon->devices;
# endif /* HW_USE_CUDA */

# ifdef HW_USE_OPENCL
    threads += opencl_recon->devices;
# endif /* HW_USE_OPENCL */

# if HW_USE_CPU	== 1
    threads += cpu_recon->devices - 1; // one core reserved for data preloading
# elif defined(HW_USE_CPU)
    if (!threads) threads += cpu_recon->devices - 1;
# endif /* HW_USE_CPU */
#else /* HW_USE_THREADS */
    threads = 1;
#endif /* HW_USE_THREADS */

    ctx->num_recon = threads;
    ctx->recon = (HSTReconstructorContext**)malloc((ctx->num_recon + 1) * sizeof(HSTReconstructorContext*));
    check_alloc(ctx->recon, "pool of HSTReconstructors");

    ctx->recon[ctx->num_recon] = NULL;

    i = 0; limit = 0;
#ifdef HW_USE_CUDA
    limit += cuda_recon->devices;
    for (; (i < limit)&&(i < threads); i++) {
	ctx->recon[i] = hst_reconstructor_create(cuda_recon, &ctx->setup, i);
	check_alloc(ctx->recon[i], "GPUContext");
    } 
#endif /* HW_USE_CUDA */

#ifdef HW_USE_OPENCL
    limit += opencl_recon->devices;
    for (; (i < limit) && (i < threads); i++) {
        ctx->recon[i] = hst_reconstructor_create(opencl_recon, &ctx->setup, i);
        check_alloc(ctx->recon[i], "OpenCLContext");
    }
#endif /* HW_USE_OPENCL */

#ifdef HW_USE_CPU
    for (; i < threads; i++) {
	ctx->recon[i] = hst_reconstructor_create(cpu_recon, &ctx->setup, threads - i - 1);
	check_alloc(ctx->recon[i], "CPUContext");
    }
#endif /* HW_USE_CPU */

    if (!i) pyhst_fail("no reconstructors found");
    
    pyhst_warning("Running %i reconstructors\n", threads);

    ctx->sched = hw_sched_create(threads);
    check_alloc(ctx->sched, "scheduller");
    
    ctx->cpu_sched = hw_sched_create((cpu_recon->devices>1)?2:1);
//    ctx->cpu_sched = hw_sched_create((cpu_recon->devices>1)?(cpu_recon->devices - 1):1);
    check_alloc(ctx->cpu_sched, "cpu scheduller");

    ctx->output_mutex = hw_sched_create_mutex();
    //Can return NULL if threads are disabled
    //check_alloc(ctx->output_mutex, "hw_sched_create_mutex");

#ifdef HST_MULTISLICE_MODE
    hw_sched_set_sequential_mode(ctx->sched, &ctx->num_slices, &ctx->current_slice, HW_SCHED_FLAG_FREE_CALL|HW_SCHED_FLAG_TERMINATOR_CALL);
#else /* HST_MULTISLICE_MODE */
    hw_sched_set_sequential_mode(ctx->sched, &ctx->num_slices, &ctx->current_slice, 0);
#endif /* HST_MULTISLICE_MODE */

    err = hw_sched_schedule_thread_task(ctx->sched, ctx, hst_threads_init);
    check_code(err, "hw_sched_schedule_task");
    
    hw_sched_set_sequential_mode(ctx->cpu_sched, &ctx->num_projections, &ctx->current_projection, 0);

/*
    for (i = 0; i < threads; i++) {
	err = hst_reconstructor_init(ctx->recon[i]);
	check_code(err, "hst_cpu_init_context");
    }
*/
    ctx->projection_data = (float**)malloc(ctx->setup.num_projections * sizeof(float*));
    check_alloc(ctx->projection_data, "projection_data");

#ifdef HST_MULTISLICE_MODE
    ctx->result_tmpbuf = hst_pinned_malloc(ctx->setup.max_slices * ctx->setup.num_x * ctx->setup.num_y * sizeof(float), HST_PINNED_W);
#else /* HST_MULTISLICE_MODE */
    ctx->result_tmpbuf = hst_pinned_malloc(threads * ctx->setup.num_x * ctx->setup.num_y * sizeof(float), HST_PINNED_W);
#endif /* HST_MULTISLICE_MODE */
    check_alloc(ctx->result_tmpbuf, "temporary buffer for reconstructed slices");

    err = hw_sched_wait_thread_task(ctx->sched);
    check_code(err, "hw_sched_wait_task");

    return 0;
}

static int hst_threads_free(HWThread thr, void *hwctx, int device_id, void *data) {
    HSTContext *ctx = (HSTContext*)data;
    hst_reconstructor_free(ctx->recon[device_id]);
    
    return 0;
}

void hst_free_context(HSTContext *ctx) {
    int i;
    
    assert(ctx);

    hst_pinned_free(ctx->result_tmpbuf);

    if (ctx->projection_data) {
        free(ctx->projection_data);
    }

    if (ctx->recon) {
	if (ctx->sched) {
	    hw_sched_execute_thread_task(ctx->sched, ctx, hst_threads_free);
	}
	
	for (i = 0; i < ctx->num_recon; i++) {
	    if (ctx->recon[i]) hst_reconstructor_destroy(ctx->recon[i]);
	}

	free(ctx->recon);
    }

    if (ctx->output_mutex) {
	hw_sched_destroy_mutex(ctx->output_mutex);
    }
    
    if (ctx->cpu_sched) {
        hw_sched_destroy(ctx->cpu_sched);
    }
    
    if (ctx->sched) {
	hw_sched_destroy(ctx->sched);
    }

    hst_setup_free(&ctx->setup);
    
    memset(ctx, 0, sizeof(HSTContext));
}

void hst_destroy_context(HSTContext *ctx) {
    --hst_copies;
    hst_free_context(ctx);
    free(ctx);
}