/tomo/pyhst

To get this branch, use:
bzr branch http://darksoft.org/webbzr/tomo/pyhst

« back to all changes in this revision

Viewing changes to hst.c

  • Committer: Suren A. Chilingaryan
  • Date: 2017-09-28 10:34:47 UTC
  • Revision ID: csa@suren.me-20170928103447-dggjgnuxmymgew2a
Quality fixes (tex-based)

Show diffs side-by-side

added added

removed removed

Lines of Context:
hst.c
217
217
        // Computing FFT dimensions, hard to understand why the -1 is there but it works nicely
218
218
        // DS: This is actually two times more than the closest power of 2. Do we need that?
219
219
    ctx->setup.dim_fft = 1<<((int)(log((1. * ctx->setup.fft_oversampling * ctx->setup.num_bins - 1)) / log(2.0) + 0.9999));
 
220
    printf("FFT convolution: %u\n", ctx->setup.dim_fft);
220
221
    ctx->setup.filter = malloc((2 * ctx->setup.dim_fft) * sizeof(float));
221
222
    check_alloc(ctx->setup.filter, " FILTER");
222
223
 
514
515
    return 0;
515
516
}
516
517
 
517
 
 
 
518
#include "astra_fourier.h"
518
519
int hst_set_filter(HSTContext *ctx, HSTFilterFunction func, void *func_param) {
519
520
    int j;
520
521
    int dim_fft;
521
522
    float *FILTER;
522
523
 
523
524
    assert(ctx);
524
 
    
 
525
 
525
526
    dim_fft = ctx->setup.dim_fft;
526
527
    FILTER = ctx->setup.filter;
527
528
 
 
529
#ifdef PYHST_ASTRA_SCALING
 
530
        // dim_fft is next power of two of (2 * num_proj)
 
531
        // num_bins is dim_fft/2 + 1
 
532
    int num_bins = dim_fft / 2 + 1;     // number of non-symmetric coefficients
 
533
 
 
534
    memset(FILTER, 0, 2 * dim_fft * sizeof(float));
 
535
    FILTER[0] = 0.25;
 
536
    FILTER[1] = 0;
 
537
    for (j = 1; j < dim_fft; j++) {
 
538
        if (j&1) {
 
539
            double coef;
 
540
            if (2 * j > dim_fft)
 
541
                coef = M_PI * (dim_fft - j);
 
542
            else
 
543
                coef = M_PI * j;
 
544
            FILTER[2 * j] = -1. / (coef * coef);
 
545
        } else {
 
546
            FILTER[2 * j] = 0;
 
547
        }
 
548
        FILTER[2 * j + 1] = 0;
 
549
    }
 
550
 
 
551
    int *ip = (int*)malloc((2 + sqrt(dim_fft) + 1)*sizeof(int)); ip[0] = 0;
 
552
    float *w = (float*)malloc(dim_fft * sizeof(float) / 2);
 
553
    cdft(2 * dim_fft, -1, FILTER, ip, w);
 
554
    free(w);
 
555
    free(ip);
 
556
 
 
557
    for (j = 0; j < num_bins; j++) {
 
558
        float _fD = 1.0f;
 
559
        float fRelIndex = (float)j / (float)dim_fft;
 
560
        float pfW = 2. * M_PI * fRelIndex;
 
561
 
 
562
            // On last iteration it is actually equal. So, rounding error works differently with ASTRA and here.... Therefore, we require a bit on top.
 
563
        if (pfW > (M_PI * _fD + 1E-6)) {
 
564
            FILTER[2 * j] = 0;
 
565
            FILTER[2 * j + 1] = 0;
 
566
        } else {
 
567
            FILTER[2 * j] = 2.0 * FILTER[2 * j];
 
568
            FILTER[2 * j + 1] = FILTER[2 * j];
 
569
        }
 
570
    }
 
571
 
 
572
    for (j = num_bins; j < dim_fft; j++) {
 
573
        FILTER[2 * j] = 0;
 
574
        FILTER[2 * j + 1] = 0;
 
575
    }
 
576
 
 
577
/*
 
578
    for (j = 0; j < num_bins; j++) {
 
579
        printf("%8.4f %8.4f ", FILTER[2 * j], FILTER[2 * j + 1]);
 
580
    }
 
581
    printf("\n");
 
582
*/
 
583
#else /* ASTRA_SCALING */
528
584
    if (func) {
529
585
        func(func_param, dim_fft, FILTER);
530
586
    } else {
531
587
        FILTER[0]=0.0;
532
 
        FILTER[1]= 1.0/dim_fft ;             /* This is the real part of the highest frequency */
 
588
        FILTER[1]= 1.0/dim_fft ;             // This is the real part of the highest frequency 
533
589
 
534
590
        for (j=1;j<dim_fft / 2; ++j) {
535
 
            FILTER[2 * j] = j * 2.0 / dim_fft/ dim_fft;
 
591
            FILTER[2 * j] = j * 2.0 / dim_fft / dim_fft;
536
592
            FILTER[2 * j + 1] = FILTER[2 * j];
537
593
        }
538
594
    }
544
600
    if (ctx->setup.no_filtering) {
545
601
        FILTER[1] = -9999.9999;
546
602
    }
547
 
 
548
603
    /* This is due to different organization of data representation used after by
549
604
    real to complex FFT transfer. In both cases only non-redundant coefficients
550
605
    are stored however, the current CPU code stores High Frequency Term (HFT) in the
551
606
    res[1], while cuFFT (and FFTW as well) have always res[1] = 0 and stores
552
607
    the HFT in res[dim_fft] */
553
608
 
554
 
    FILTER[dim_fft] = FILTER[1];
 
609
    FILTER[dim_fft] = FILTER[1]; // or should it be 0?
555
610
 
556
611
    /* CUDA code now includes optimization allowing to compute two convolutions 
557
612
    parallely using single Complex-Complex Fourier transform. For that we need
562
617
 
563
618
    FILTER[dim_fft + 1] = FILTER[1];
564
619
    FILTER[1] = FILTER[0];
565
 
 
 
620
#endif /* ASTRA_SCALING */
566
621
 
567
622
    for (j = dim_fft + 2; j < 2 * dim_fft; j+=2) {
568
623
        FILTER[j] = FILTER[2 * dim_fft - j];