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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
|
/*
-----------------------------------------------------------------------
Copyright 2012 iMinds-Vision Lab, University of Antwerp
Contact: astra@ua.ac.be
Website: http://astra.ua.ac.be
This file is part of the
All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox").
The ASTRA Toolbox 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.
The ASTRA Toolbox 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 the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
-----------------------------------------------------------------------
$Id$
*/
#ifndef _CUDA_ASTRA_H
#define _CUDA_ASTRA_H
#include "fft.h"
#include "fbp_filters.h"
#include "dims.h"
#include "algo.h"
using astraCUDA::SFanProjection;
namespace astra {
enum Cuda2DProjectionKernel {
ker2d_default = 0
};
class AstraFBP_internal;
class _AstraExport AstraFBP {
public:
// Constructor
AstraFBP();
// Destructor
~AstraFBP();
// Set the size of the reconstruction rectangle.
// Volume pixels are currently assumed to be 1x1 squares.
bool setReconstructionGeometry(unsigned int iVolWidth,
unsigned int iVolHeight,
float fPixelSize = 1.0f);
// Set the projection angles and number of detector pixels per angle.
// pfAngles must be a float array of length iProjAngles.
// fDetSize indicates the size of a detector pixel compared to a
// volume pixel edge.
//
// pfAngles will only be read from during this call.
bool setProjectionGeometry(unsigned int iProjAngles,
unsigned int iProjDets,
const float *pfAngles,
float fDetSize = 1.0f);
// Set linear supersampling factor for the BP.
// (The number of rays is the square of this)
//
// This may optionally be called before init().
bool setPixelSuperSampling(unsigned int iPixelSuperSampling);
// Set per-detector shifts.
//
// pfTOffsets will only be read from during this call.
bool setTOffsets(const float *pfTOffsets);
// Returns the required size of a filter in the fourier domain
// when multiplying it with the fft of the projection data.
// Its value is equal to the smallest power of two larger than
// or equal to twice the number of detectors in the spatial domain.
//
// _iDetectorCount is the number of detectors in the spatial domain.
static int calcFourierFilterSize(int _iDetectorCount);
// Sets the filter type. Some filter types require the user to supply an
// array containing the filter.
// The number of elements in a filter in the fourier domain should be equal
// to the value returned by calcFourierFilterSize().
// The following types require a filter:
//
// - FILTER_PROJECTION:
// The filter size should be equal to the output of
// calcFourierFilterSize(). The filtered sinogram is
// multiplied with the supplied filter.
//
// - FILTER_SINOGRAM:
// Same as FILTER_PROJECTION, but now the filter should contain a row for
// every projection direction.
//
// - FILTER_RPROJECTION:
// The filter should now contain one kernel (= ifft of filter), with the
// peak in the center. The filter width
// can be any value. If odd, the peak is assumed to be in the center, if
// even, it is assumed to be at floor(filter-width/2).
//
// - FILTER_RSINOGRAM
// Same as FILTER_RPROJECTION, but now the supplied filter should contain a
// row for every projection direction.
//
// A large number of other filters (FILTER_RAMLAK, FILTER_SHEPPLOGAN,
// FILTER_COSINE, FILTER_HAMMING, and FILTER_HANN)
// have a D variable, which gives the cutoff point in the frequency domain.
// Setting this value to 1.0 will include the whole filter
bool setFilter(E_FBPFILTER _eFilter,
const float * _pfHostFilter = NULL,
int _iFilterWidth = 0, float _fD = 1.0f, float _fFilterParameter = -1.0f);
// Initialize CUDA, allocate GPU buffers and
// precompute geometry-specific data.
//
// CUDA is set up to use GPU number iGPUIndex.
//
// This must be called after calling setReconstructionGeometry() and
// setProjectionGeometry().
bool init(int iGPUIndex = 0);
// Setup input sinogram for a slice.
// pfSinogram must be a float array of size iProjAngles*iSinogramPitch.
// NB: iSinogramPitch is measured in floats, not in bytes.
//
// This must be called after init(), and before iterate(). It may be
// called again after iterate()/getReconstruction() to start a new slice.
//
// pfSinogram will only be read from during this call.
bool setSinogram(const float* pfSinogram, unsigned int iSinogramPitch);
// Runs an FBP reconstruction.
// This must be called after setSinogram().
//
// run can be called before setFilter, but will then use the default Ram-Lak filter
bool run();
// Get the reconstructed slice.
// pfReconstruction must be a float array of size
// iVolHeight*iReconstructionPitch.
// NB: iReconstructionPitch is measured in floats, not in bytes.
//
// This may be called after run().
bool getReconstruction(float* pfReconstruction,
unsigned int iReconstructionPitch) const;
private:
AstraFBP_internal* pData;
};
class _AstraExport BPalgo : public astraCUDA::ReconAlgo {
public:
BPalgo();
~BPalgo();
virtual bool init();
virtual bool iterate(unsigned int iterations);
virtual float computeDiffNorm();
};
// TODO: Clean up this interface to FP
// Do a single forward projection
_AstraExport bool astraCudaFP(const float* pfVolume, float* pfSinogram,
unsigned int iVolWidth, unsigned int iVolHeight,
unsigned int iProjAngles, unsigned int iProjDets,
const float *pfAngles, const float *pfOffsets,
float fDetSize = 1.0f, unsigned int iDetSuperSampling = 1,
int iGPUIndex = 0);
// Do a single forward projection, fan beam
_AstraExport bool astraCudaFanFP(const float* pfVolume, float* pfSinogram,
unsigned int iVolWidth, unsigned int iVolHeight,
unsigned int iProjAngles, unsigned int iProjDets,
const float *pfAngles, float fOriginSourceDistance,
float fOriginDetectorDistance, float fPixelSize = 1.0f,
float fDetSize = 1.0f,
unsigned int iDetSuperSampling = 1,
int iGPUIndex = 0);
_AstraExport bool astraCudaFanFP(const float* pfVolume, float* pfSinogram,
unsigned int iVolWidth, unsigned int iVolHeight,
unsigned int iProjAngles, unsigned int iProjDets,
const SFanProjection *pAngles,
unsigned int iDetSuperSampling = 1,
int iGPUIndex = 0);
}
#endif
|
