/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);