1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
|
/*
* 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 <math.h>
#include "hst_reconstructor.h"
int hst_reconstructor_init_context(HSTReconstructorContext *ctx, HSTReconstructor *recon, HSTSetup *setup) {
memcpy(ctx, recon, sizeof(HSTReconstructor));
ctx->setup = setup;
ctx->processed_slices = 0;
return 0;
}
void hst_reconstructor_free_context(HSTReconstructorContext *ctx) {
}
int hst_reconstructor_postprocess_slice(HSTReconstructorContext *ctx, float *slice, const float *sinograms) {
HSTSetup *setup;
int i, x, y;
int num_x;
int num_y;
int num_bins;
int num_projections;
/* The block of variables to reproduce the fundamental, needed if SUMRULE = 1*/
float Sino_Sum = 0;
float Slice_Sum = 0;
assert(ctx);
setup = ctx->setup;
assert(setup);
num_x = setup->num_x;
num_y = setup->num_y;
num_bins = setup->num_bins;
num_projections = setup->num_projections;
if (setup->zerooffmask) {
assert(setup->minX);
assert(setup->maxX);
for (y = 0; y < num_y; y++) {
for (x = 0; x <= setup->minX[y]; x++) {
slice[y*num_x +x] = 0;
}
for (x = setup->maxX[y]; x < num_x; x++) {
slice[y*num_x +x] = 0;
}
}
}
if (setup->sum_rule == 1) {
/* This could be moved to GPU (see reduction sample in SDK). */
for (i = 0; i < num_x * num_y; i++) {
Slice_Sum += slice[i];
}
for (i = 0; i < num_projections * num_bins;i++) {
Sino_Sum += sinograms[i];
}
Sino_Sum *= 2.0 / M_PI ;
for (i = 0; i < num_x*num_y; i++) {
slice[i] += (Sino_Sum - Slice_Sum) / (num_x * num_y);
}
}
return 0;
}
|