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
|
#include "filter_convolve.h"
#include "complexdata.h"
#include "fitting.h"
class FilterFindFWHM : public MinimizationFunction {
public:
FilterFindFWHM(const LDRfilter& kernel) : kern(kernel) {}
virtual unsigned int numof_fitpars() const {return 1;}
virtual float evaluate(const fvector& fv) const {return fabs(kern.calculate(fv[0])-0.5);}
private:
const LDRfilter& kern;
};
////////////////////////////////////////////////////////////////////////////////
void FilterConvolve::init() {
kernel.set_description("convolution kernel");
append_arg(kernel,"kernel");
kernelwidth.set_unit(ODIN_SPAT_UNIT).set_description("full-width-at-half-maximum of kernel");
append_arg(kernelwidth,"kernelwidth");
}
bool FilterConvolve::process(Data<float,4>& data, Protocol& prot) const {
Log<Filter> odinlog(c_label(),"process");
Range all=Range::all();
FilterFindFWHM findfwhm(kernel);
DownhillSimplex downhill(findfwhm);
fvector starting_point(1); starting_point[0]=0.5;
fvector step_size(1); step_size[0]=0.1;
fvector fitresult(1);
float relfwhm=1.0;
if(downhill.get_minimum_parameters(fitresult, starting_point, step_size)) {
relfwhm=fitresult[0];
}
float kerneldiameter=secureDivision(kernelwidth,relfwhm);
ODINLOG(odinlog,normalDebug) << "kernelwidth/relfwhm/kerneldiameter=" << kernelwidth << "/" << relfwhm << "/" << kerneldiameter << STD_endl;
TinyVector<int,4> shape=data.shape();
TinyVector<int,3> spatshape(shape(sliceDim), shape(phaseDim), shape(readDim));
TinyVector<float,3> voxel_spacing;
voxel_spacing(0)=FileFormat::voxel_extent(prot.geometry, sliceDirection, data.extent(sliceDim));
voxel_spacing(1)=FileFormat::voxel_extent(prot.geometry, phaseDirection, data.extent(phaseDim));
voxel_spacing(2)=FileFormat::voxel_extent(prot.geometry, readDirection, data.extent(readDim));
ODINLOG(odinlog,normalDebug) << "voxel_spacing=" << voxel_spacing << STD_endl;
int oversampling=2;
TinyVector<int,3> kernelshape=oversampling*spatshape;
TinyVector<float,3> voxel_spacing_kernelgrid(voxel_spacing/oversampling);
ODINLOG(odinlog,normalDebug) << "voxel_spacing_kernelgrid=" << voxel_spacing_kernelgrid << STD_endl;
ComplexData<3> kernelgrid(kernelshape);
kernelgrid=STD_complex(0.0);
// fill with gridding kernel
for(unsigned int i=0; i<kernelgrid.size(); i++) {
TinyVector<int,3> index=kernelgrid.create_index(i);
TinyVector<float,3> dist=voxel_spacing_kernelgrid*(index-kernelshape/2);
float relradius=secureDivision( sqrt(double(sum(dist*dist))), 0.5*kerneldiameter);
if(relradius<=1.0) {
float kernelval=kernel.calculate(relradius);
kernelgrid(index)=STD_complex(kernelval);
}
}
kernelgrid.fft(true);
// Data<float,3>(cabs(kernelgrid)).autowrite("kernelgrid.nii");
TinyVector<int,3> lowerBounds(kernelshape/4);
TinyVector<int,3> upperBounds(lowerBounds+spatshape-1);
ComplexData<3> kernelmult(spatshape);
kernelmult=kernelgrid(RectDomain<3>(lowerBounds, upperBounds));
float kernelmax=max(creal(kernelmult));
ODINLOG(odinlog,normalDebug) << "kernelmax=" << kernelmax << STD_endl;
if(kernelmax>0.0) kernelmult=kernelmult/STD_complex(kernelmax);
// Data<float,3>(cabs(kernelmult)).autowrite("kernelmult.nii");
for(int irep=0; irep<shape(timeDim); irep++) {
ComplexData<3> onevol(float2real(data(irep,all,all,all)));
onevol.fft(true);
onevol=onevol*kernelmult;
onevol.fft(false);
data(irep,all,all,all)=creal(onevol);
}
return true;
}
|
