gridding.h
1 /***************************************************************************
2  gridding.h - description
3  -------------------
4  begin : Fri Jun 25 2004
5  copyright : (C) 2004 by Michael von Mengershausen
6  email : mengers@cns.mpg.de
7  ***************************************************************************/
8 
9 /***************************************************************************
10  * *
11  * This program is free software; you can redistribute it and/or modify *
12  * it under the terms of the GNU General Public License as published by *
13  * the Free Software Foundation; either version 2 of the License, or *
14  * (at your option) any later version. *
15  * *
16  ***************************************************************************/
17 
18 #ifndef GRIDDING_H
19 #define GRIDDING_H
20 
21 #include <odindata/data.h>
22 #include <odindata/utils.h>
23 
24 #include <odinpara/ldrfilter.h>
25 
36 template<int N_rank>
37 struct GriddingPoint {
38  GriddingPoint(const TinyVector<float,N_rank>& c=0.0, float w=1.0) : coord(c), weight(w) {}
39 
40  TinyVector<float,N_rank> coord;
41  float weight;
42 };
43 
44 
46 
51 template<typename T, int N_rank>
52 class Gridding {
53 
54  public:
55 
59  Gridding() : shape(0) {}
60 
70  Array<float,N_rank> init(const TinyVector<int,N_rank>& dst_shape, const TinyVector<float,N_rank>& dst_extent,
71  const STD_vector<GriddingPoint<N_rank> >& src_coords,
72  const LDRfilter& kernel, float kernel_diameter);
73 
80  template<int N_rank_in>
81  Array<T,N_rank> operator () (const Array<T,N_rank_in>& src, unsigned int offset=0) const;
82 
83 
84  private:
85  TinyVector<int,N_rank> shape;
86 
87  // numof src-steps x numof indices on dst grid x pair(dst-index, weight)
88  STD_vector< STD_vector< STD_pair<TinyVector<int,N_rank>, float> > > recipe;
89 
90 
91 };
92 
94 
95 template <typename T, int N_rank>
96 Array<float,N_rank> Gridding<T,N_rank>::init(const TinyVector<int,N_rank>& dst_shape, const TinyVector<float,N_rank>& dst_extent,
97  const STD_vector<GriddingPoint<N_rank> >& src_coords,
98  const LDRfilter& kernel, float kernel_diameter) {
99  Log<OdinData> odinlog("Gridding","init");
100 
101  shape=dst_shape;
102 
103  unsigned int nsrc=src_coords.size();
104 
105  ODINLOG(odinlog,normalDebug) << "nsrc/kernel/kernel_diameter=" << nsrc << "/" << kernel.get_function_name() << "/" << kernel_diameter << STD_endl;
106 
107  recipe.resize(nsrc);
108 
109  // array to collect sum of all weights for later normalization
110  Array<float,N_rank> weight_sum(dst_shape);
111  weight_sum=0.0;
112 
113  TinyVector<float,N_rank> dst_step=dst_extent/dst_shape;
114  ODINLOG(odinlog,normalDebug) << "dst_step=" << dst_step << STD_endl;
115 
116  TinyVector<float,N_rank> kernel_extent; // In units of indices
117  for(int irank=0; irank<N_rank; irank++) {
118  if(dst_step(irank)>0.0) kernel_extent(irank)=kernel_diameter/dst_step(irank);
119  else kernel_extent(irank)=0.0;
120  }
121  ODINLOG(odinlog,normalDebug) << "kernel_extent=" << kernel_extent << STD_endl;
122 
123  TinyVector<float,N_rank> offset=0.5*(dst_shape-1.0); // coordinate is at center of destination voxels
124 
125  // Iterate over src coordinates
126  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
127  const GriddingPoint<N_rank>& point=src_coords[isrc];
128 
129  // Find grid points in neighboordhood
130  TinyVector<float,N_rank> root; // In units of indices on destination grid
131  for(int irank=0; irank<N_rank; irank++) {
132  if(dst_step(irank)>0.0) root(irank)=point.coord(irank)/dst_step(irank);
133  else root(irank)=0.0;
134  }
135  root += offset;
136  ODINLOG(odinlog,normalDebug) << "root=" << root << STD_endl;
137 
138  TinyVector<int,N_rank> lowindex=(root-0.5*kernel_extent)+0.5; // round up to next higher grid point
139  TinyVector<int,N_rank> uppindex=(root+0.5*kernel_extent); // round down to next lower grid point
140  ODINLOG(odinlog,normalDebug) << "lowindex/uppindex=" << lowindex << "/" << uppindex << STD_endl;
141 
142  // Possible points in neighboorhood
143  TinyVector<int,N_rank> neighbours=uppindex-lowindex+1;
144  ODINLOG(odinlog,normalDebug) << "neighbours=" << neighbours << STD_endl;
145 
146  // Get actual points in neighbourhood and their weight if they are within the support of the given kernel and whether point lies on destination grid
147  STD_vector<STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[isrc];
148  dstvec.clear();
149  for(int ineighb=0; ineighb<product(neighbours); ineighb++) {
150  TinyVector<int,N_rank> neighb_index=index2extent(neighbours, ineighb);
151 
152  // Check whether point is on grid
153  TinyVector<int,N_rank> index=lowindex+neighb_index;
154  bool ongrid=true;
155  for(int i=0; i<N_rank; i++) {
156  if(index(i)<0) ongrid=false;
157  if(index(i)>=dst_shape(i)) ongrid=false;
158  }
159 
160  // Check whether point is within the support of the kernel
161  if(ongrid) {
162  TinyVector<float,N_rank> diff=(root-index)*dst_step;
163  float radius=sqrt(sum(diff*diff))/(0.5*kernel_diameter);
164  float weight=point.weight*kernel.calculate(radius);
165  if(weight>=0.0) dstvec.push_back(STD_pair<TinyVector<int,N_rank>, float>(index,weight));
166  }
167  }
168 
169  // Sum up weights for each index
170  for(unsigned int idst=0; idst<dstvec.size(); idst++) {
171  weight_sum(dstvec[idst].first)+=dstvec[idst].second;
172  }
173 
174  } // End of iterating over src coordinates
175 
176 
177  // Normalize weights
178  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
179  STD_vector< STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[isrc];
180  for(unsigned int idst=0; idst<dstvec.size(); idst++) {
181  STD_pair<TinyVector<int,N_rank>, float>& weightpair=dstvec[idst];
182  float weightsum=weight_sum(weightpair.first);
183  if(weightsum>0.0) weightpair.second/=weightsum;
184  }
185  }
186 
187  return weight_sum;
188 }
189 
191 
192 template <typename T, int N_rank>
193 template<int N_rank_in>
194 Array<T,N_rank> Gridding<T,N_rank>::operator () (const Array<T,N_rank_in>& src, unsigned int offset) const {
195  Log<OdinData> odinlog("Gridding","()");
196  Array<T,N_rank> result;
197 
198  unsigned int nsrc=src.size();
199  if((offset+nsrc)>recipe.size()) {
200  ODINLOG(odinlog,errorLog) << "Max index of src=" << offset+nsrc << " exceeds recipe.size()=" << recipe.size() << STD_endl;
201  return result;
202  }
203 
204  result.resize(shape);
205  result=T(0);
206 
207  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
208  const STD_vector< STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[offset+isrc];
209  for(unsigned int idst=0; idst<dstvec.size(); idst++) {
210  const STD_pair<TinyVector<int,N_rank>, float>& weightpair=dstvec[idst];
211  result(weightpair.first) += weightpair.second * src(index2extent(src.shape(),isrc));
212  }
213  }
214 
215 
216  return result;
217 }
218 
219 
220 
222 
226 template<typename T, int N_rank, bool OnPixelRot=false>
228 
229  public:
230 
236  CoordTransformation(const TinyVector<int,N_rank>& shape, const TinyMatrix<float,N_rank,N_rank>& rotation, const TinyVector<float,N_rank>& offset, float kernel_size=2.5);
237 
241  Array<T,N_rank> operator () (const Array<T,N_rank>& A) const;
242 
243  private:
244  TinyVector<int,N_rank> shape_cache;
245  Gridding<T,N_rank> gridder;
246 };
247 
249 
250 template <typename T, int N_rank, bool OnPixelRot>
251 CoordTransformation<T,N_rank,OnPixelRot>::CoordTransformation(const TinyVector<int,N_rank>& shape, const TinyMatrix<float,N_rank,N_rank>& rotation, const TinyVector<float,N_rank>& offset, float kernel_size)
252  : shape_cache(shape) {
253  Log<OdinData> odinlog("CoordTransformation","CoordTransformation");
254 
255  int nsrc=product(shape);
256  STD_vector<GriddingPoint<N_rank> > src_coords(nsrc);
257  ODINLOG(odinlog,normalDebug) << "N_rank/nsrc=" << N_rank << "/" << nsrc << STD_endl;
258 
259  TinyVector<float,N_rank> extent_half;
260  if(OnPixelRot) extent_half=shape/2; // suitable for rotating k-space
261  else extent_half=0.5*(shape-1); // center of even size is between voxels
262  TinyVector<int,N_rank> index;
263  TinyVector<float,N_rank> findex;
264  for(int isrc=0; isrc<nsrc; isrc++) {
265  index=index2extent(shape, isrc);
266 
267  // calculating new coordinates
268  findex=index-extent_half;
269  src_coords[isrc].coord=rotation*findex+offset; // Modify coordinate of source points, this will effectively rotate/shift the image data after gridding
270  }
271  ODINLOG(odinlog,normalDebug) << "source done" << STD_endl;
272 
273 
274  LDRfilter gridkernel;
275  gridkernel.set_function("Gauss");
276  gridder.init(shape, shape, src_coords, gridkernel, kernel_size);
277 }
278 
280 
281 
282 template <typename T, int N_rank, bool OnPixelRot>
283 Array<T,N_rank> CoordTransformation<T,N_rank,OnPixelRot>::operator () (const Array<T,N_rank>& A) const {
284  Log<OdinData> odinlog("CoordTransformation","()");
285  if(A.shape()!=shape_cache) {
286  ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
287  return A;
288  }
289 
290  return gridder(A);
291 }
292 
293 
299 #endif
Gridding()
Definition: gridding.h:59
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
Definition: tjlog.h:218
Array< T, N_rank > operator()(const Array< T, N_rank_in > &src, unsigned int offset=0) const
Definition: gridding.h:194
Array< T, N_rank > operator()(const Array< T, N_rank > &A) const
Definition: gridding.h:283
CoordTransformation(const TinyVector< int, N_rank > &shape, const TinyMatrix< float, N_rank, N_rank > &rotation, const TinyVector< float, N_rank > &offset, float kernel_size=2.5)
Definition: gridding.h:251