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
|
#include "filter_align.h"
#include "gridding.h"
void FilterAlign::init(){
fname.set_description("filename");
append_arg(fname,"fname");
blowup.set_description("In-plane blowup factor");
append_arg(blowup,"blowup");
}
bool FilterAlign::process(Data<float,4>& data, Protocol& prot) const {
Log<Filter> odinlog(c_label(),"process");
Range all=Range::all();
TinyVector<int,3> inshape3d(data.extent(sliceDim),data.extent(phaseDim),data.extent(readDim));
int blowup_factor=STD_max(1,int(blowup));
// Load external file
Data<float,4> extdata;
Protocol alignprot;
if(extdata.autoread(fname, FileReadOpts(), &alignprot)<0) return false;
TinyVector<int,3> dst_shape(extdata.shape()(sliceDim),blowup_factor*extdata.shape()(phaseDim),blowup_factor*extdata.shape()(readDim));
ODINLOG(odinlog,normalDebug) << "dst_shape=" << dst_shape << STD_endl;
extdata.free();
TinyVector<float,3> dst_extent;
for(int i=0; i<3; i++) dst_extent(2-i)=alignprot.geometry.get_FOV(direction(i));
ODINLOG(odinlog,normalDebug) << "dst_extent=" << dst_extent << STD_endl;
int npts3d=product(inshape3d);
ODINLOG(odinlog,normalDebug) << "npts3d=" << npts3d << STD_endl;
STD_vector<GriddingPoint<3> > src_coords(npts3d);
dvector srcfov(3);
dvector dstfov(3);
dvector srcsize(3);
dvector dstsize(3);
for(int i=0; i<n_directions; i++) {
srcfov[i]=prot.geometry.get_FOV(direction(i));
dstfov[i]=alignprot.geometry.get_FOV(direction(i));
srcsize[i]=inshape3d(2-i);
dstsize[i]=dst_shape(2-i);
}
ODINLOG(odinlog,normalDebug) << "srcfov/dstfov=" << srcfov << "/" << dstfov << STD_endl;
ODINLOG(odinlog,normalDebug) << "srcsize/dstsize=" << srcsize << "/" << dstsize << STD_endl;
// Create src coordinates
dvector relpos(3);
dvector srcpos(3);
dvector dstpos(3);
for(int ipt=0; ipt<npts3d; ipt++) {
TinyVector<int,3> index=index2extent(inshape3d, ipt);
for(int i=0; i<3; i++) relpos[2-i]=secureDivision(0.5+index(i), inshape3d(i));
srcpos=srcfov*(relpos-0.5);
dstpos=alignprot.geometry.transform( prot.geometry.transform(srcpos,false), true);
src_coords[ipt].coord=TinyVector<float,3>(dstpos[2],dstpos[1],dstpos[0]);
ODINLOG(odinlog,normalDebug) << "src_coords(" << ipt << ")=" << src_coords[ipt].coord << STD_endl;
}
LDRfilter gridkernel;
gridkernel.set_function("Gauss");
// Calculate kernel which covers both grids
for(int i=0; i<n_directions; i++) {
if(srcsize[i]<2) srcfov[i]=0.0;
if(dstsize[i]<2) dstfov[i]=0.0;
}
dvector srcspacing(srcfov/srcsize);
dvector dstspacing(dstfov/dstsize);
float kernel_diameter=STD_max( sqrt((srcspacing*srcspacing).sum()), sqrt((dstspacing*dstspacing).sum()) );
ODINLOG(odinlog,normalDebug) << "kernel_diameter=" << kernel_diameter << STD_endl;
Gridding<float,3> gridder;
gridder.init(dst_shape, dst_extent, src_coords, gridkernel, kernel_diameter);
int ntime=data.extent(timeDim);
Data<float,4> outdata(ntime,dst_shape(0),dst_shape(1),dst_shape(2));
for(int itime=0; itime<data.extent(timeDim); itime++) {
#ifdef BLITZ_REPLACEMENT
outdata(itime,all,all,all)=gridder.operator()<3>(data(itime,all,all,all)); // syntax for blitz_replacement
#else
outdata(itime,all,all,all)=gridder(data(itime,all,all,all));
#endif
}
data.reference(outdata);
// Copy relevant part of protocol
prot.geometry=alignprot.geometry;
for(int i=0; i<n_directions; i++) prot.seqpars.set_MatrixSize(direction(i),alignprot.seqpars.get_MatrixSize(direction(i)));
return true;
}
|
