/tomo/pyhst

/*
 * The PyHST program is Copyright (C) 2002-2008 of the
 * European Synchrotron Radiation Facility (ESRF).
 *
 * 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 <math.h>

#include <glib.h>

#ifdef PYHST_USE_FFTW
# include <fftw3.h>
#else /* PYHST_USE_FFTW */
# include "Vhst_fourier.h"
#endif /* PYHST_USE_FFTW */

#include "hw_sched.h"
#include "Vhst_calculate_limits.h"

#include "cpumain.h"

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

#include "hst_reconstructor.h"



static HSTConstString hst_cpu_timers[] = { "complete", "initialization", "filtering", "backprojection", NULL };

#define CPU_CONTEXT(ctx) ((CPUContext*)ctx)

struct CPUContextT {
    HSTReconstructorContext recon;
	
#ifdef PYHST_MEASURE_TIMINGS
# define CPU_CONTEXT_MEMSET_OFFSET (3 * sizeof(GTimer*))
    GTimer *main_timer;		// Counts time spent in backprojection code
    GTimer *pre_timer;		// Counts time spent in filtering
    GTimer *init_timer;		// Counts time spent in memory operations
#else
# define CPU_CONTEXT_MEMSET_OFFSET 0
#endif /* PYHST_MEASURE_TIMINGS */

    float *WORK;			// Preallocated temporary buffer
    float *OVERSAMPLE;		// Bigger temporary buffer for interpolation

#ifdef PYHST_USE_FFTW
    HWSched sched;
    fftwf_plan forward_plan;
    fftwf_plan inverse_plan;
#else /* PYHST_USE_FFTW */
    int dim_exponts;		// dimension of EXPONENTS array
    float *EXPONENTS;		// EXPONENTS used for FFT computation by Vhst_fourier implementation
#endif /* PYHST_USE_FFTW */
};
typedef struct CPUContextT CPUContext;


#ifndef PYHST_USE_FFTW
static int hst_cpu_compute_exponents(int dim_exponts,int dim_fft2,int radix,float *EXPONENTS) {
    int status = 0;
    hst_expnt__(&dim_exponts,&dim_fft2,&radix,EXPONENTS,&status);
    return status;
}
#endif /* !PYHST_USE_FFTW */

static HSTReconstructorContext *hst_cpu_create_context(HSTReconstructor *prototype, HSTSetup *setup, int id) {
    assert(setup);
    
    CPUContext *ctx = (CPUContext*)malloc(sizeof(CPUContext));

    if (ctx) {
        memset(ctx, 0, sizeof(CPUContext));

#ifdef PYHST_MEASURE_TIMINGS
        ctx->main_timer = g_timer_new();
        if (ctx->main_timer) g_timer_stop(ctx->main_timer);
        ctx->pre_timer = g_timer_new();
        if (ctx->pre_timer) g_timer_stop(ctx->pre_timer);
        ctx->init_timer = g_timer_new();
        if (ctx->init_timer) g_timer_stop(ctx->init_timer);

        if ((!ctx->main_timer)||(!ctx->pre_timer)||(!ctx->init_timer)) {
            if (ctx->main_timer) g_timer_destroy(ctx->main_timer);
            if (ctx->pre_timer) g_timer_destroy(ctx->pre_timer);
            if (ctx->init_timer) g_timer_destroy(ctx->init_timer);
            free(ctx);
            return NULL;
        }
#endif /* PYHST_MEASURE_TIMINGS */

	hst_reconstructor_init_context((HSTReconstructorContext*)ctx, prototype, setup);

    }
    return (HSTReconstructorContext*)ctx;
}

static int hst_cpu_init_context(HSTReconstructorContext *rctx, HWThread thr) {
    HSTSetup *setup;
    CPUContext *ctx;

    assert(rctx);
    

    ctx = CPU_CONTEXT(rctx);
    setup = rctx->setup;

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_continue(ctx->init_timer);
#endif /* PYHST_MEASURE_TIMINGS */
    
    ctx->WORK = malloc((setup->dim_fft + 2) * sizeof(float));
    check_alloc(ctx->WORK,"WORK" );
    memset(ctx->WORK,0, (setup->dim_fft + 2) * sizeof(float));

    /*******************************************************************************
     * Takes care of oversampling space
     ******************************************************************************/
    ctx->OVERSAMPLE = (float*)malloc(setup->num_projections * 3 * setup->num_bins * setup->oversampling * sizeof(float)) ;
    check_alloc(ctx->OVERSAMPLE,"OVERSAMPLE");
    memset(ctx->OVERSAMPLE , 0, 3*setup->num_projections*setup->num_bins*setup->oversampling*sizeof(float));

    /*******************************************************************************
     * Takes care of fourier space
     ******************************************************************************/
#ifdef PYHST_USE_FFTW
    ctx->sched = thr->sched;
    hw_sched_lock(ctx->sched, sync);
    printf("dim_fft %lu\n", setup->dim_fft);
    ctx->forward_plan = fftwf_plan_dft_r2c_1d(setup->dim_fft, ctx->WORK, (fftwf_complex*)ctx->WORK, 0);
    ctx->inverse_plan = fftwf_plan_dft_c2r_1d(setup->dim_fft, (fftwf_complex*)ctx->WORK, ctx->WORK, 0);
    hw_sched_unlock(ctx->sched, sync);
#else /* PYHST_USE_FFTW */
    ctx->dim_exponts = 2 * setup->dim_fft;
    ctx->EXPONENTS = malloc(ctx->dim_exponts * sizeof(float));
    check_alloc(ctx->EXPONENTS,"EXPONENTS");

    /************************************************************
     * calculation of exponential factors for fft              *
     ************************************************************/
    hst_cpu_compute_exponents(ctx->dim_exponts, setup->dim_fft*2, 8, ctx->EXPONENTS);
#endif /* PYHST_USE_FFTW */

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_stop(ctx->init_timer);
#endif /* PYHST_MEASURE_TIMINGS */

    return 0;
}

static void hst_cpu_free_context(HSTReconstructorContext *rctx) {
    CPUContext *ctx;

    ctx = CPU_CONTEXT(rctx);

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_continue(ctx->init_timer);
#endif /* PYHST_MEASURE_TIMINGS */

#ifdef PYHST_USE_FFTW
    if (ctx->sched) {
	hw_sched_lock(ctx->sched, sync);
	if (ctx->inverse_plan) fftwf_destroy_plan(ctx->inverse_plan);
	if (ctx->forward_plan) fftwf_destroy_plan(ctx->forward_plan);
	hw_sched_unlock(ctx->sched, sync);
    }
#else /* PYHST_USE_FFTW */
    if (ctx->EXPONENTS) free(ctx->EXPONENTS);
#endif /* PYHST_USE_FFTW */
    if (ctx->WORK) free(ctx->WORK);
    if (ctx->OVERSAMPLE) free(ctx->OVERSAMPLE);
    memset(((char*)ctx) + sizeof(HSTReconstructorContext) + CPU_CONTEXT_MEMSET_OFFSET, 0, sizeof(CPUContext) - sizeof(HSTReconstructorContext) - CPU_CONTEXT_MEMSET_OFFSET);

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_stop(ctx->init_timer);
#endif /* PYHST_MEASURE_TIMINGS */

}

static void hst_cpu_destroy_context(HSTReconstructorContext *rctx) {
#ifdef PYHST_MEASURE_TIMINGS
    CPUContext *ctx;
#endif /* PYHST_MEASURE_TIMINGS */

    hst_cpu_free_context(rctx);
    hst_reconstructor_free_context(rctx);

#ifdef PYHST_MEASURE_TIMINGS
    ctx = CPU_CONTEXT(rctx);

    g_timer_destroy(ctx->pre_timer);
    g_timer_destroy(ctx->main_timer);
    g_timer_destroy(ctx->init_timer);
#endif /* PYHST_MEASURE_TIMINGS */

    free(rctx);
}

static HSTConstString hst_cpu_get_title(HSTReconstructorConstContextPtr rctx) {
    return "CPU";
}

static HSTConstString *hst_cpu_get_timers(HSTReconstructorConstContextPtr rctx, double *timers) {
#ifdef PYHST_MEASURE_TIMINGS
    assert(rctx);
    
    if (timers) {
	timers[1] = g_timer_elapsed(CPU_CONTEXT(rctx)->init_timer, NULL);
	timers[2] = g_timer_elapsed(CPU_CONTEXT(rctx)->pre_timer, NULL);
	timers[3] = g_timer_elapsed(CPU_CONTEXT(rctx)->main_timer, NULL);
	timers[0] = timers[1] + timers[2] + timers[3];
    }
#endif /* PYHST_MEASURE_TIMINGS */

    return hst_cpu_timers;
}


static int hst_cpu_preprocess(HSTReconstructorContext *rctx, float *SLICE, const float *SINOGRAMS) {
    // The block of general iterators
    int i,j;
    long projection;

#ifndef PYHST_USE_FFTW
    // Temporary variables
    int status = 0;
    int True = 1;
#endif /* !PYHST_USE_FFTW */

    // Variables needed for applying limit in fai360 mode
    float dx;
    int prof_shift = 0, prof_fact = 0;
    float overlapping = 0, flat_zone = 0;
    float pente_zone;

    // Copies
    CPUContext *ctx;
    HSTSetup *setup;
    float *WORK;
    float *OVERSAMPLE;
    float *FILTER;
    int dim_fft, num_bins, num_projections;
    float axis_position, axis_position_corr;
    float *axis_position_corr_s;
    int oversampling;

    ctx = CPU_CONTEXT(rctx);

    assert(ctx);
    assert(ctx->WORK);
    assert(ctx->OVERSAMPLE);

#ifdef PYHST_USE_FFTW
    assert(ctx->forward_plan);
    assert(ctx->inverse_plan);
#else /* PYHST_USE_FFTW */
    assert(ctx->EXPONENTS);
#endif /* PYHST_USE_FFTW */

    setup = rctx->setup;
    WORK = ctx->WORK;
    OVERSAMPLE = ctx->OVERSAMPLE;
    FILTER = setup->filter;
    num_projections = setup->num_projections;
    num_bins = setup->num_bins;
    dim_fft = setup->dim_fft;
    axis_position = setup->axis_position;
    axis_position_corr_s = setup->axis_position_corr_s;
    oversampling = setup->oversampling;

    pente_zone = setup->pente_zone;


    for (projection = 0; projection < num_projections; projection++) {
        axis_position_corr = axis_position_corr_s[projection];
        memcpy(WORK, SINOGRAMS + projection * num_bins, num_bins * sizeof(float));

    	// PADDING
        if (setup->zero_padding) {
            memset(WORK + num_bins, 0, (dim_fft - num_bins) * sizeof(float));
        } else {
            for (i = num_bins ; i < (num_bins + dim_fft) / 2 ; i++) {
                WORK[i] = WORK[num_bins -1];
            }
            for (i = (num_bins + dim_fft) / 2 ; i < dim_fft ; i++) {
                WORK[i] = WORK[0];
            }
        }

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_continue(ctx->pre_timer);
#endif /* PYHST_MEASURE_TIMINGS */

#ifdef PYHST_USE_FFTW
	fftwf_execute(ctx->forward_plan);
#else /* PYHST_USE_FFTW */
        hst_ffa8__(&dim_fft, &dim_fft, &ctx->dim_exponts, ctx->EXPONENTS, &True, WORK, &status);
#endif /* PYHST_USE_FFTW */
        /*
        {
        FILE *f = fopen("/tmp/fourier.txt", "a+");
        for (i = 0; i < dim_fft + 2; i++) fprintf(f, "%6u - %+.8e,", i, WORK[i]);
        fprintf(f,"\n");
        fclose(f);
        }
        */

        /******************************************************
        * Multiply the Fourier transformed array with filter *
        ******************************************************/
if (0) {
#ifdef PYHST_USE_FFTW
        for (i = 0;i <= dim_fft;i++) {
#else /* PYHST_USE_FFTW */
	    // Embeded FFT engine have slightly different data organization, see hst_set_filter for explanations
	WORK[0] = WORK[0] * FILTER[0];
	WORK[1] = WORK[1] * FILTER[dim_fft];
        for (i = 2;i < dim_fft;i++) {
#endif /* PYHST_USE_FFTW */
            WORK[i] = WORK[i] * FILTER[i]; // fft e inv lasciano un fattore dim_fft, filter nella forma a rampa arriva fino a 1/dim_fft 
        }
}
        /***************************************************
        * Inverse Fourier transform to give convolution *
        ***************************************************/

#ifdef PYHST_USE_FFTW
	fftwf_execute(ctx->inverse_plan);
#else /* PYHST_USE_FFTW */
        hst_ffs8__(&dim_fft, &dim_fft, &ctx->dim_exponts, ctx->EXPONENTS, &True, WORK, &status);
#endif /* PYHST_USE_FFTW */

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_stop(ctx->pre_timer);
#endif /* PYHST_MEASURE_TIMINGS */


        if (setup->fai360) {
            if (2*axis_position_corr > num_bins) {
                overlapping = (num_bins - axis_position_corr);
                prof_fact = 1;
                prof_shift = 0;

            } else {
                overlapping = (axis_position_corr);
                prof_fact = -1;
                prof_shift = 1;
            }

            if (overlapping <= 0) {
                pyhst_fail("corrected axis falls out of CCD");
            }

            pente_zone = MIN(pente_zone, overlapping);
            flat_zone = overlapping - pente_zone;

            for (i = 0; i < num_bins;i++) {
                if (fabs(i - axis_position_corr) <= overlapping) {
                    if (fabs(i - axis_position_corr) <= flat_zone) {
                        WORK[i] *= 0.5;
                    } else {
                        if (i > axis_position_corr) {
                            dx = axis_position_corr - i + flat_zone;
                        } else {
                            dx = axis_position_corr - i - flat_zone;
                        }
                        WORK[i] *= prof_shift + prof_fact * (0.5 * (1.0 + (dx) / pente_zone));
                    }
                } else if (i > axis_position_corr) {
                    WORK[i] *= prof_shift - 0.0;
                } else {
                    WORK[i] *= prof_shift + prof_fact * 1.0 ;
                }
                WORK[i] *= 2;
            }
        }

        {
            WORK[num_bins] = WORK[num_bins-1];
            for (i = 0; i < num_bins; i++) {
                for (j = 0; j < oversampling; j++) {
                    OVERSAMPLE[ (projection*3 + 1)*num_bins*oversampling + i*oversampling + j ] =
                        (WORK[i] * (oversampling - j) + WORK[i+1] * j) / oversampling;
                }
            }
        }
    } // end of projections loop


    return 0;
}


static int hst_cpu_reconstruct(HSTReconstructorContext *rctx, float *SLICE, const float *SINOGRAMS) {
    CPUContext *ctx;
    HSTSetup *setup;

    ctx = CPU_CONTEXT(rctx);

    assert(ctx);

    setup = rctx->setup;

    memset(SLICE, 0, setup->num_x * setup->num_y * sizeof(float));


/*
    Just for testing: original version of PyHST have a bug by not always setting axis_position_corr_s (for
    test example it is only 900), so just set it here to get completely indentical input for cpu_main.
    setup->axis_position_corr_s[900]=0;
*/
/*
    FILE *f;
    f = fopen("/tmp/orig_redesig.txt", "w");
    fwrite (ctx->OVERSAMPLE,4,3*(setup->num_projections)*(setup->num_bins)*setup->oversampling, f);
    fwrite(setup->axis_position_corr_s,4,setup->num_projections,f);
    fwrite(setup->cos_s,4,setup->num_projections,f);
    fwrite(setup->sin_s,4,setup->num_projections,f);
    fwrite(SLICE,4,setup->num_x*setup->num_y,f);
    fclose(f);
    printf("%u %u %u %u %u\n", setup->num_y, setup->num_x, setup->num_projections, setup->num_bins, setup->oversampling);
    printf("%f %f %f\n", setup->axis_position, setup->offset_x, setup->offset_y);
*/

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_continue(ctx->main_timer);
#endif /* PYHST_MEASURE_TIMINGS */

    cpu_main(setup->num_y, setup->num_x, SLICE, setup->num_projections, setup->num_bins, ctx->OVERSAMPLE,
             setup->axis_position, setup->axis_position_corr_s, setup->cos_s, setup->sin_s,
             setup->offset_x, setup->offset_y, NULL, NULL, /*minX, maxX,*/ setup->oversampling);

#ifdef PYHST_MEASURE_TIMINGS
    g_timer_stop(ctx->main_timer);
#endif /* PYHST_MEASURE_TIMINGS */

    return 0;
}


static HSTReconstructor hst_cpu_info = {
    1,
    hst_cpu_get_title,
    hst_cpu_create_context,
    hst_cpu_init_context,
    hst_cpu_free_context,
    hst_cpu_destroy_context,
    NULL,
    NULL,
    hst_cpu_preprocess,
    hst_cpu_reconstruct,
    hst_reconstructor_postprocess_slice,
    NULL,
    hst_cpu_get_timers
};


HSTReconstructor *hst_cpu_init(int flags) {
    // Detect number of cores (embedded FFT is thread-unsafe)
#ifdef PYHST_USE_FFTW
    hst_cpu_info.devices = hw_sched_get_cpu_count()/2;
#endif /* PYHST_USE_FFTW */

    return &hst_cpu_info;
}

void hst_cpu_free() {
}