gridding.h

/***************************************************************************
                          gridding.h  -  description
                             -------------------
    begin                : Fri Jun 25 2004
    copyright            : (C) 2004-2021 by Michael von Mengershausen and Thies Jochimsen
    email                : mengers@cns.mpg.de
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program 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 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#ifndef GRIDDING_H
#define GRIDDING_H
 
#include <odindata/data.h>
#include <odindata/utils.h>
 
#include <odinpara/ldrfilter.h>
 
template<int N_rank>
struct GriddingPoint {
  GriddingPoint(const TinyVector<float,N_rank>& c=0.0, float w=1.0) : coord(c), weight(w) {}
 
  TinyVector<float,N_rank> coord;
  float weight;
};
 
 
 
template<typename T, int N_rank>
class Gridding {
 
 public:
 
  Gridding() : shape(0) {}
 
  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);
 
  template<int N_rank_in>
  Array<T,N_rank> operator () (const Array<T,N_rank_in>& src, unsigned int offset=0) const;
 
 
 private:
  TinyVector<int,N_rank> shape;
 
  // numof src-steps  x  numof indices on dst grid  x  pair(dst-index, weight)
  STD_vector< STD_vector< STD_pair<TinyVector<int,N_rank>, float> > > recipe;
 
 
};
 
 
template <typename T, int N_rank>
Array<float,N_rank> Gridding<T,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) {
 Log<OdinData> odinlog("Gridding","init");
 
  shape=dst_shape;
 
  unsigned int nsrc=src_coords.size();
 
  ODINLOG(odinlog,normalDebug) << "nsrc/kernel/kernel_diameter=" << nsrc << "/" << kernel.get_function_name() << "/" << kernel_diameter << STD_endl;
 
  recipe.resize(nsrc);
 
  // array to collect sum of all weights for later normalization
  Array<float,N_rank> weight_sum(dst_shape);
  weight_sum=0.0;
 
  TinyVector<float,N_rank> dst_step=dst_extent/dst_shape;
  ODINLOG(odinlog,normalDebug) << "dst_step=" << dst_step << STD_endl;
 
  TinyVector<float,N_rank> kernel_extent; // In units of indices
  for(int irank=0; irank<N_rank; irank++) {
    if(dst_step(irank)>0.0) kernel_extent(irank)=kernel_diameter/dst_step(irank);
    else kernel_extent(irank)=0.0;
  }
  ODINLOG(odinlog,normalDebug) << "kernel_extent=" << kernel_extent << STD_endl;
 
  TinyVector<float,N_rank> offset=0.5*(dst_shape-1.0); // coordinate is at center of destination voxels
 
  // Iterate over src coordinates
  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
    const GriddingPoint<N_rank>& point=src_coords[isrc];
 
    // Find grid points in neighboordhood
    TinyVector<float,N_rank> root; // In units of indices on destination grid
    for(int irank=0; irank<N_rank; irank++) {
      if(dst_step(irank)>0.0) root(irank)=point.coord(irank)/dst_step(irank);
      else root(irank)=0.0;
    }
    root += offset;
    ODINLOG(odinlog,normalDebug) << "root=" << root << STD_endl;
 
    TinyVector<int,N_rank> lowindex=(root-0.5*kernel_extent)+0.5; // round up to next higher grid point
    TinyVector<int,N_rank> uppindex=(root+0.5*kernel_extent); // round down to next lower grid point
    ODINLOG(odinlog,normalDebug) << "lowindex/uppindex=" << lowindex << "/" << uppindex << STD_endl;
 
    // Possible points in neighboorhood
    TinyVector<int,N_rank> neighbours=uppindex-lowindex+1;
    ODINLOG(odinlog,normalDebug) << "neighbours=" << neighbours << STD_endl;
 
    // 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
    STD_vector<STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[isrc];
    dstvec.clear();
    for(int ineighb=0; ineighb<product(neighbours); ineighb++) {
      TinyVector<int,N_rank> neighb_index=index2extent(neighbours, ineighb);
 
      // Check whether point is on grid
      TinyVector<int,N_rank> index=lowindex+neighb_index;
      bool ongrid=true;
      for(int i=0; i<N_rank; i++) {
        if(index(i)<0)             ongrid=false;
        if(index(i)>=dst_shape(i)) ongrid=false;
      }
 
      // Check whether point is within the support of the kernel
      if(ongrid) {
        TinyVector<float,N_rank> diff=(root-index)*dst_step;
        float radius=sqrt(sum(diff*diff))/(0.5*kernel_diameter);
        float weight=point.weight*kernel.calculate(radius);
        if(weight>=0.0) dstvec.push_back(STD_pair<TinyVector<int,N_rank>, float>(index,weight));
      }
    }
 
    // Sum up weights for each index
    for(unsigned int idst=0; idst<dstvec.size(); idst++) {
      weight_sum(dstvec[idst].first)+=dstvec[idst].second;
    }
 
  } // End of iterating over src coordinates
 
 
  // Normalize weights
  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
    STD_vector< STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[isrc];
    for(unsigned int idst=0; idst<dstvec.size(); idst++) {
      STD_pair<TinyVector<int,N_rank>, float>& weightpair=dstvec[idst];
      float weightsum=weight_sum(weightpair.first);
      if(weightsum>0.0) weightpair.second/=weightsum;
    }
  }
 
  return weight_sum;
}
 
 
template <typename T, int N_rank>
template<int N_rank_in>
Array<T,N_rank> Gridding<T,N_rank>::operator () (const Array<T,N_rank_in>& src, unsigned int offset) const {
  Log<OdinData> odinlog("Gridding","()");
  Array<T,N_rank> result;
 
  unsigned int nsrc=src.size();
  if((offset+nsrc)>recipe.size()) {
    ODINLOG(odinlog,errorLog) << "Max index of src=" << offset+nsrc << " exceeds recipe.size()=" << recipe.size() << STD_endl;
    return result;
  }
 
  result.resize(shape);
  result=T(0);
 
  for(unsigned int isrc=0; isrc<nsrc; isrc++) {
    const STD_vector< STD_pair<TinyVector<int,N_rank>, float> >& dstvec=recipe[offset+isrc];
    for(unsigned int idst=0; idst<dstvec.size(); idst++) {
      const STD_pair<TinyVector<int,N_rank>, float>& weightpair=dstvec[idst];
      result(weightpair.first) += weightpair.second * src(index2extent(src.shape(),isrc));
    }
  }
 
 
  return result;
}
 
 
 
 
template<typename T, int N_rank, bool OnPixelRot=false>
class CoordTransformation {
 
 public:
 
  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);
 
  Array<T,N_rank> operator () (const Array<T,N_rank>& A) const;
 
 private:
  TinyVector<int,N_rank> shape_cache;
  Gridding<T,N_rank> gridder;
};
 
 
template <typename T, int N_rank, bool OnPixelRot>
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)
 : shape_cache(shape) {
  Log<OdinData> odinlog("CoordTransformation","CoordTransformation");
 
  int nsrc=product(shape);
  STD_vector<GriddingPoint<N_rank> > src_coords(nsrc);
  ODINLOG(odinlog,normalDebug) << "N_rank/nsrc=" << N_rank << "/" << nsrc << STD_endl;
 
  TinyVector<float,N_rank> extent_half;
  if(OnPixelRot) extent_half=shape/2; // suitable for rotating k-space
  else extent_half=0.5*(shape-1);     // center of even size is between voxels
  TinyVector<int,N_rank> index;
  TinyVector<float,N_rank> findex;
  for(int isrc=0; isrc<nsrc; isrc++) {
    index=index2extent(shape, isrc);
 
    // calculating new coordinates
    findex=index-extent_half;
    src_coords[isrc].coord=rotation*findex+offset;  // Modify coordinate of source points, this will effectively rotate/shift the image data after gridding
  }
  ODINLOG(odinlog,normalDebug) << "source done" << STD_endl;
 
 
  LDRfilter gridkernel;
  gridkernel.set_function("Gauss");
  gridder.init(shape, shape, src_coords, gridkernel, kernel_size);
}
 
 
 
template <typename T, int N_rank, bool OnPixelRot>
Array<T,N_rank> CoordTransformation<T,N_rank,OnPixelRot>::operator () (const Array<T,N_rank>& A) const {
  Log<OdinData> odinlog("CoordTransformation","()");
  if(A.shape()!=shape_cache) {
    ODINLOG(odinlog,errorLog) << "Shape mismatch" << STD_endl;
    return A;
  }
 
  return gridder(A);
}
 
 
#endif