grid_code.h
1 #include "grid.h"
2 #include "data.h"
3 #include "controller.h"
4 
5 
6 static const char* deapodize_post_str="deapodize";
7 
8 
9 template<int NDim>
10 void RecoGrid<NDim>::init() {
11 
12  deapo_shift=true;
13  deapo_shift.set_description("Shift deapodization function for correct offset");
14  append_arg(deapo_shift,"deapo_shift");
15 
16  kernel.set_funcpars("Gauss");
17  kernel.set_cmdline_option("kg"+postfix()).set_description("Kernel function for "+postfix()+" k-space gridding");
18  append_arg(kernel,"kernel");
19 
20  kernelwidth=sqrt(float(3-NDim))*3.0; // found empirically
21  kernelwidth.set_cmdline_option("kw"+postfix()).set_unit("Pixel").set_description("Kernel width for "+postfix()+" k-space gridding");
22  append_arg(kernelwidth,"kernelwidth");
23 
24  gridScale = 1.0;
25  gridScale.set_cmdline_option("gs").set_description("Scale Factor for the grid used for gridding");
26  append_arg(gridScale,"gridScale");
27 }
28 
29 
30 
31 template<int NDim>
32 TinyVector<int,NDim> RecoGrid<NDim>::gridshape(const RecoCoord& coord, RecoController& controller) {
33  Log<Reco> odinlog("RecoGrid","gridshape");
34 
35  TinyVector<int,3> imagesize=controller.image_size();
36  TinyVector<int,NDim> result;
37  for(int i=0; i<NDim; i++) result(i)=imagesize(3-NDim+i);
38 
39  if(NDim==1) {
40  if(coord.readoutShape && (coord.readoutShape->second)>0) result(0)=coord.readoutShape->second; // already includes oversampling
41  else result(0)=int(result(0)*coord.overSamplingFactor+0.5); // only for ADC regridding, oversampling will be removed later in image domain
42  }
43  ODINLOG(odinlog,normalDebug) << "result=" << result << STD_endl;
44 
45  return result;
46 }
47 
48 
49 
50 template<int NDim>
52  Log<Reco> odinlog(c_label(),"get_gridder");
53 
54  MutexLock lock(gridmutex);
55 
56  const KspaceTraj* traj=coord.kspaceTraj; // access is protected by mutex
57 
58  typename GridderMap::const_iterator it=gridmap.find(traj);
59  if(it==gridmap.end()) {
60  Gridding<STD_complex,NDim>& gridder=gridmap[traj];
61 
62  float appGridScale = gridScale;
63  if (NDim == 1) appGridScale = 1.0; // don't use grid scaling for 1D Gridding
64 
65  TinyVector<int,NDim> dstshape=gridshape(coord, controller);
66 
67  int nseg=traj->extent(0); // e.g. spiral segments
68  int npts=traj->extent(1);
69 
70  STD_vector<GriddingPoint<NDim> > srccoords(nseg*npts); // all segments in one gridder for correct weighting of src values
71 
72  for(int iseg=0; iseg<nseg; iseg++) {
73  for(int ipt=0; ipt<npts; ipt++) {
74  const GriddingPoint<3>& gridpoint=(*traj)(iseg,ipt);
75  const TinyVector<float,3>& kpoint=gridpoint.coord;
76  GriddingPoint<NDim>& srccoord=srccoords[iseg*npts+ipt];
77  TinyVector<float,NDim>& srcpoint=srccoord.coord;
78  for(int i=0; i<NDim; i++) {
79  float k=kpoint(3-NDim+i); // Use rightmost k-space coordinate subset
80  srcpoint(i)=k * appGridScale;
81  }
82  ODINLOG(odinlog,normalDebug) << "srcpoint=" << srcpoint << STD_endl;
83  ODINLOG(odinlog,normalDebug) << "kpoint=" << kpoint << STD_endl;
84  srccoord.weight=gridpoint.weight;
85  }
86  }
87 
88  TinyVector<float,3> kmax=controller.kspace_extent();
89  ODINLOG(odinlog,normalDebug) << "kmax=" << kmax << STD_endl;
90 
91  TinyVector<float,NDim> dstextent;
92  for(int i=0; i<NDim; i++) dstextent(i)=kmax(3-NDim+i);
93 
94  for (int i=0; i<NDim; i++) {
95  dstshape(i) *= appGridScale;
96  dstextent(i) *= appGridScale;
97  }
98 
99  float dk=max(dstextent/dstshape);
100  float kwidth=kernelwidth*dk*appGridScale;
101  ODINLOG(odinlog,normalDebug) << "dstshape/dstextent/dk/kwidth=" << dstshape << "/" << dstextent << "/" << dk << "/" << kwidth << STD_endl;
102 
103  Data<float,NDim> sampldens(dstshape);
104  sampldens=gridder.init(dstshape, dstextent, srccoords, kernel, kwidth);
105 // sampldens.autowrite("sampldens.jdx");
106 
107 
108 
109  // Create magnitude correction for regridded data according to the gridding kernel used
110  int deapo_os_factor=4; // oversampling to account for small kernels
111  TinyVector<int,NDim> shape_os(deapo_os_factor*dstshape);
112  ComplexData<NDim> kspace_kernel_os(shape_os);
113  kspace_kernel_os=STD_complex(0.0);
114 
115  // fill k-space with gridding kernel
116  for(unsigned int i=0; i<kspace_kernel_os.size(); i++) {
117  TinyVector<int,NDim> index=kspace_kernel_os.create_index(i);
118  TinyVector<int,NDim> dist=index-shape_os/2;
119  float radius=secureDivision( sqrt(double(sum(dist*dist))), 0.5*deapo_os_factor*kernelwidth*appGridScale);
120  if(radius<=1.0) kspace_kernel_os(index)=STD_complex(kernel.calculate(radius));
121  }
122 
123 
124  ODINLOG(odinlog,normalDebug) << "deapo_shift=" << deapo_shift << STD_endl;
125 
126  if(deapo_shift) {
127  // Take spatial offset into account
128  TinyVector<float,3> reloffset=controller.reloffset();
129  TinyVector<float,NDim> reloffset_os;
130  for(int i=0; i<NDim; i++) reloffset_os(i)=reloffset(3-NDim+i)/deapo_os_factor/appGridScale;
131  if(max(abs(reloffset))>0.0) kspace_kernel_os.modulate_offset(reloffset_os);
132  }
133 
134  // FFT and cutting out central part with de-apodize function
135  kspace_kernel_os.fft();
136  Data<float,NDim> apodize(dstshape);
137  apodize=cabs(kspace_kernel_os(RectDomain<NDim>((shape_os-dstshape)/2, (shape_os-dstshape)/2+dstshape-1)));
138 // apodize.autowrite("apodize"+ptos(traj)+".jdx");
139 
140  if(min(apodize)>0.0) {
141  RecoData rdposted;
142 
143  // post separately for templates
144  rdposted.coord().index[templtype]=coord.index[templtype];
145  rdposted.coord().index[templtype].set_mode(RecoIndex::separate);
146 
147  rdposted.data(Rank<NDim>()).reference(float2real(mean(apodize)/apodize));
148  controller.post_data(deapodize_post_str+postfix(), rdposted);
149  } else {
150  ODINLOG(odinlog,warningLog) << "FT of gridding kernel contains zero, not applying de-apodize correction" << STD_endl;
151  }
152 
153 
154  return gridder;
155  }
156  return it->second;
157 }
158 
159 
160 
161 template<int NDim>
162 bool RecoGrid<NDim>::process(RecoData& rd, RecoController& controller) {
163  Log<Reco> odinlog(c_label(),"process");
164 
165  if(!rd.coord().kspaceTraj) {
166  ODINLOG(odinlog,errorLog) << "No kspaceTraj available for coord=" << rd.coord().print() << STD_endl;
167  return false;
168  }
169 
170  rd.coord().gridScaleFactor = gridScale;
171 
172  const Gridding<STD_complex,NDim>& gridder=get_gridder(rd.coord(), controller);
173 
174  TinyVector<int,NDim> outshape=gridshape(rd.coord(), controller);
175  if (NDim > 1) {
176  for (int i=0; i<NDim; i++) outshape(i) *= gridScale;
177  }
178 
179  ComplexData<NDim> outdata(outshape);
180 
181 
182  ComplexData<1>& adc=rd.data(Rank<1>());
183  int adcSize=adc.extent(0);
184  ODINLOG(odinlog,normalDebug) << "adcSize=" << adcSize << STD_endl;
185 
186  unsigned int offset=0;
187  if(rd.coord().kspaceTraj->extent(0)>1) offset=rd.coord().index[cycle]*adcSize; // only for segmented 2D spirals, ignore cycle otherwise
188 
189  outdata=gridder(adc, offset);
190 
191  rd.data(Rank<NDim>()).reference(outdata);
192 
193  return execute_next_step(rd,controller);
194 }
195 
196 
197 template<int NDim>
199  if(context.mode==RecoQueryContext::prep) {
200  context.controller.announce_data(deapodize_post_str+postfix());
201  }
202  return RecoStep::query(context);
203 }
204 
205 
207 
208 template<int NDim>
210 
211  RecoData rddeapo;
212  if(!controller.inquire_data(deapodize_post_str+postfix(), rddeapo)) return false;
213 
214  ComplexData<3>& data=rd.data(Rank<3>());
215  ComplexData<NDim>& deapo=rddeapo.data(Rank<NDim>());
216 
217  for(unsigned int i=0; i<data.size(); i++) {
218  TinyVector<int,3> index=data.create_index(i);
219 
220  TinyVector<int,NDim> deapoindex;
221  for(int i=0; i<NDim; i++) deapoindex(i)=index(3-NDim+i);
222 
223  data(index)*=deapo(deapoindex);
224 
225  }
226 
227 
228 // for(int iphase3d=0; iphase3d<data.extent(0); iphase3d++) data(iphase3d,all,all)=data(iphase3d,all,all)*deapo(all,all);
229 
230  return execute_next_step(rd,controller);
231 }
232 
233 
234 
236 
237 bool RecoGridCut::process(RecoData& rd, RecoController& controller) {
238 
239  float scaleFactor = rd.coord().gridScaleFactor;
240  if (scaleFactor != 1.0) {
241 
242  Range all=Range::all();
243 
244  ComplexData<3>& indata=rd.data(Rank<3>());
245  TinyVector<int,3> inshape(indata.shape());
246  TinyVector<int,3> outshape(inshape(0), inshape(1) / scaleFactor, inshape(2) / scaleFactor);
247 
248  TinyVector<int,2> phasebounds((inshape(1) - outshape(1)) / 2, (inshape(1) + outshape(1)) / 2);
249  TinyVector<int,2> readbounds((inshape(1) - outshape(1)) / 2, (inshape(1) + outshape(1)) / 2);
250 
251  ComplexData<3> outdata(outshape);
252  outdata = indata(all,Range(phasebounds(0),phasebounds(1)),Range(readbounds(0),readbounds(1)));
253  rd.data(Rank<3>()).reference(outdata);
254 
255  rd.coord().gridScaleFactor = 1.0; // reset to 1.0 since the scaling has been removed from the data
256  }
257 
258  return execute_next_step(rd,controller);
259 }
260 
261 
TinyVector< int, N_rank > create_index(unsigned long index) const
bool execute_next_step(RecoData &rd, RecoController &controller)
virtual bool query(RecoQueryContext &context)
TinyVector< int, 3 > image_size() const
STD_string print(RecoIndex::printMode printmode=RecoIndex::brief) const
Array< float, N_rank > init(const TinyVector< int, N_rank > &dst_shape, const TinyVector< float, N_rank > &dst_extent, const STD_vector< GriddingPoint< N_rank > > &src_coords, const LDRfilter &kernel, float kernel_diameter)
Definition: gridding.h:96
const KspaceTraj * kspaceTraj
Definition: index.h:289
Definition: tjlog.h:218
RecoCoord & coord()
RecoIndex & set_mode(indexMode m)
Definition: index.h:80
bool inquire_data(const STD_string &label, RecoData &data)
TinyVector< float, 3 > reloffset() const
TinyVector< float, 3 > kspace_extent() const
const ReadoutShape * readoutShape
Definition: index.h:282
grid1d, grid2d: Grid signal data with arbitrary k-space trajectory in 2 dimensions ...
Definition: grid.h:26
RecoIndex index[n_recoDims]
Definition: index.h:269
void announce_data(const STD_string &label)
Definition: controller.h:109
void post_data(const STD_string &label, const RecoData &data)
Definition: controller.h:114
float gridScaleFactor
Definition: index.h:325
deapodize1d, deapodize2d: De-apodize gridded data in 2 dimensions
Definition: grid.h:61
float overSamplingFactor
Definition: index.h:305
RecoController & controller
Definition: odinreco/step.h:55
double secureDivision(double numerator, double denominator)
ComplexData< 1 > & data(Rank< 1 >) const
Definition: odinreco/data.h:91