/perf/fdk

#define _POSIX_C_SOURCE 200112L

#define _XOPEN_SOURCE 600

#include <stdio.h>

#include <stdlib.h>

#include <malloc.h>

#include <errno.h>

#include <sys/time.h>

#include <ctype.h>

#include <stdarg.h>

#include <math.h>

#include <pthread.h>

#include <stdint.h>

#include <string.h>

#include <limits.h>

#include <getopt.h>

#include "process.h"

extern int counter;

enum {

    OPT_SLICE = 's',

    OPT_HELP

};

static struct option long_options[] = {

    {"help",			no_argument,		0, OPT_HELP },

    {"slice",			required_argument, 	0, OPT_SLICE },

    { 0, 0, 0, 0 }

};

static void Usage(int argc, char *argv[], const char *format, ...) {

    if (format) {

        va_list ap;

        va_start(ap, format);

        printf("Error %i: ", errno);

        vprintf(format, ap);

        printf("\n");

        va_end(ap);

        printf("\n");

    }

    printf(

"Usage:\n"

" %s [options]\n"

"\n"

"Options:\n"

" -s [slice]		- slice number\n"

"\n\n",

argv[0]);

    exit(0);

}

static int isnumber(const char *str) {

    int i = 0;

    for (i = 0; str[i]; i++)

        if (!isdigit(str[i])) return 0;

    return 1;

}

void main(int argc, char *argv[]) {

    int n_threads = 1;

    int alignment = 256;

    long n_proj = 3600; 

    long n_elements = 2000;

    float cor = 1060. / 16;

    float d = 1060;

    float detector_size = 400;

    float pixel_size = 400. / 2000;

    float cor_offset = 0;

    float detector_offset_u = 0, detector_offset_v = 0, detector_offset_z = 0;

    float source[]  = {0, 0, 0};

    float u_detector[] = {1, 0, 0};

    float v_detector[] = {0, 1, 0};

    float n_detector[] = {0, 0, 1};

    float *out_volume, *projections, *slice_x, *slice_y, *slice_coord_z;

    int err;

    int i, j;

    int tnum, m;

    struct thread_info *t_info;

    void *res;

    struct timeval start, end;

    long duration;

    err = posix_memalign((void**)&out_volume, alignment, n_elements * n_elements * n_elements * sizeof(float));

    if (!err) err = posix_memalign((void**)&projections, alignment, n_elements * n_elements * n_proj * sizeof(float));

    if (!err) err = posix_memalign((void**)&slice_x, alignment, n_elements * n_elements * sizeof(float));

    if (!err) err = posix_memalign((void**)&slice_y, alignment, n_elements * n_elements * sizeof(float));

    if (!err) err = posix_memalign((void**)&slice_coord_z, alignment, n_elements * sizeof(float));

    if (err) {

        printf("Can not allocate memory");

	exit(-1);

    }

    for (i = 0; i < n_elements; i++) {

	slice_coord_z[i] = (pixel_size / 16) * ((i + 1) - (n_elements + 1.) / 2);

	for (j = 0; j < n_elements; j++) {

	    slice_x[i * n_elements + j] = (pixel_size / 16) * ((i + 1) - (n_elements + 1.)/2);

	    slice_y[i * n_elements + j] = (pixel_size / 16) * ((n_elements - j) - (n_elements + 1.)/2);

	}

    }

    /* allocate memory for threads */

    t_info = (struct thread_info*) calloc(n_threads, sizeof(struct thread_info));

    if (t_info == NULL)

    {

        printf("Structure t_info memory allocation error\n");

        exit(-1);

    }

        /* fill tread_info structure and create threads */

    counter = n_elements / 2;

    unsigned char c;

    while ((c = getopt_long(argc, argv, "s:", long_options, NULL)) != (unsigned char)-1) {

	extern int optind;

	switch (c) {

         case OPT_HELP:

            Usage(argc, argv, NULL);

            break;

	 case OPT_SLICE:

            if ((!isnumber(optarg))||(sscanf(optarg, "%i", &counter) != 1)||(counter < 0)||(counter >= n_elements)) {

                Usage(argc, argv, "Invalid slice is specified (%s)", optarg);

            }

	    break;

	 default:

	    Usage(argc, argv,  "Unknown option (%s) with argument (%s)", optarg?argv[optind-2]:argv[optind-1], optarg?optarg:"(null)");

	}

    }

    gettimeofday(&start, NULL);

    for (tnum = 0; tnum < n_threads; tnum++)

    {

	t_info[tnum].thread_num = tnum;                      /* Application-defined thread # */

        t_info[tnum].n_elements = n_elements;

        t_info[tnum].n_proj = n_proj;

        t_info[tnum].cor = cor;

        t_info[tnum].d = d;

        t_info[tnum].detector_size = detector_size; 

        t_info[tnum].pixel_size = pixel_size;

        t_info[tnum].cor_offset = cor_offset;

        t_info[tnum].detector_offset_u = detector_offset_u;

        t_info[tnum].detector_offset_v = detector_offset_v;

        t_info[tnum].detector_offset_z = detector_offset_z;

        t_info[tnum].source = source;

        t_info[tnum].u_detector = u_detector; 

        t_info[tnum].v_detector = v_detector; 

        t_info[tnum].n_detector = n_detector;

        t_info[tnum].projections = projections;

        t_info[tnum].slice_x = slice_x;

        t_info[tnum].slice_y = slice_y; 

        t_info[tnum].slice_coord_z = slice_coord_z;

        t_info[tnum].out_volume = out_volume;

        t_info[tnum].slices_per_iter = n_elements;

	if (n_threads > 1) {

	    m = pthread_create(&t_info[tnum].thread_id, NULL, &process, &t_info[tnum]);

	    if (m != 0)

	    {

		printf("Cannot create thread with error %i\n", m);

        	exit(-1);

	    }

	} else {

	    process(&t_info[tnum]);

	}

    }

    if (n_threads > 1) {

	/* join threads */

	for (tnum = 0; tnum < n_threads; tnum++)

	{

	    m = pthread_join(t_info[tnum].thread_id, &res);

	    if (m != 0)

    	    {

        	printf("Cannot join threads with error %i\n", m);

        	exit(-1);

    	    }

    	}

    }

    gettimeofday(&end, NULL);

    duration = (end.tv_sec - start.tv_sec) * 1000000 + end.tv_usec - start.tv_usec;

    printf("Time: %2u:%2u:%2u (%6u)\n", duration/3600000000, (duration%3600000000)/60000000, (duration%60000000)/1000000, (duration%1000000)); 

    free(t_info);

    free(slice_coord_z);

    free(slice_y);

    free(slice_x);

    free(projections);

    free(out_volume);

}