gridding.h

/***************************************************************************
                          gridding.h  -  description
                             -------------------
    begin                : Fri Jun 25 2004
    copyright            : (C) 2004 by Michael von Mengershausen
    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/jdxfilter.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 JDXfilter& 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 JDXfilter& 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=kernel_diameter/dst_step; // In units of indices
  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=point.coord/dst_step + offset; // In units of indices on destination grid
    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=product(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;


  JDXfilter 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

---