grid_code.h

#include "grid.h"
#include "data.h"
#include "controller.h"
 
 
static const char* deapodize_post_str="deapodize";
 
 
template<int NDim>
void RecoGrid<NDim>::init() {
 
  deapo_shift=true;
  deapo_shift.set_description("Shift deapodization function for correct offset");
  append_arg(deapo_shift,"deapo_shift");
 
  kernel.set_funcpars("Gauss");
  kernel.set_cmdline_option("kg"+postfix()).set_description("Kernel function for "+postfix()+" k-space gridding");
  append_arg(kernel,"kernel");
 
  kernelwidth=sqrt(float(3-NDim))*3.0; // found empirically
  kernelwidth.set_cmdline_option("kw"+postfix()).set_unit("Pixel").set_description("Kernel width for "+postfix()+" k-space gridding");
  append_arg(kernelwidth,"kernelwidth");
 
  gridScale = 1.0;
  gridScale.set_cmdline_option("gs").set_description("Scale Factor for the grid used for gridding");
  append_arg(gridScale,"gridScale");
}
 
 
 
template<int NDim>
TinyVector<int,NDim> RecoGrid<NDim>::gridshape(const RecoCoord& coord, RecoController& controller) {
  Log<Reco> odinlog("RecoGrid","gridshape");
 
  TinyVector<int,3> imagesize=controller.image_size();
  TinyVector<int,NDim> result;
  for(int i=0; i<NDim; i++) result(i)=imagesize(3-NDim+i);
 
  if(NDim==1) {
    if(coord.readoutShape && (coord.readoutShape->second)>0) result(0)=coord.readoutShape->second; // already includes oversampling
    else result(0)=int(result(0)*coord.overSamplingFactor+0.5); // only for ADC regridding, oversampling will be removed later in image domain
  }
  ODINLOG(odinlog,normalDebug) << "result=" << result << STD_endl;
 
  return result;
}
 
 
 
template<int NDim>
const Gridding<STD_complex,NDim>& RecoGrid<NDim>::get_gridder(const RecoCoord& coord, RecoController& controller) {
  Log<Reco> odinlog(c_label(),"get_gridder");
 
  MutexLock lock(gridmutex);
 
  const KspaceTraj* traj=coord.kspaceTraj; // access is protected by mutex
 
  typename GridderMap::const_iterator it=gridmap.find(traj);
  if(it==gridmap.end()) {
    Gridding<STD_complex,NDim>& gridder=gridmap[traj];
 
    float appGridScale = gridScale;
    if (NDim == 1) appGridScale = 1.0; // don't use grid scaling for 1D Gridding
 
    TinyVector<int,NDim> dstshape=gridshape(coord, controller);
 
    int nseg=traj->extent(0); // e.g. spiral segments
    int npts=traj->extent(1);
 
    STD_vector<GriddingPoint<NDim> > srccoords(nseg*npts); // all segments in one gridder for correct weighting of src values
 
    for(int iseg=0; iseg<nseg; iseg++) {
      for(int ipt=0; ipt<npts; ipt++) {
        const GriddingPoint<3>& gridpoint=(*traj)(iseg,ipt);
        const TinyVector<float,3>& kpoint=gridpoint.coord;
        GriddingPoint<NDim>& srccoord=srccoords[iseg*npts+ipt];
        TinyVector<float,NDim>& srcpoint=srccoord.coord;
        for(int i=0; i<NDim; i++) {
          float k=kpoint(3-NDim+i); // Use rightmost k-space coordinate subset
          srcpoint(i)=k * appGridScale;
        }
        ODINLOG(odinlog,normalDebug) << "srcpoint=" << srcpoint << STD_endl;
        ODINLOG(odinlog,normalDebug) << "kpoint=" << kpoint << STD_endl;
        srccoord.weight=gridpoint.weight;
      }
    }
 
    TinyVector<float,3> kmax=controller.kspace_extent();
    ODINLOG(odinlog,normalDebug) << "kmax=" << kmax << STD_endl;
 
    TinyVector<float,NDim> dstextent;
    for(int i=0; i<NDim; i++) dstextent(i)=kmax(3-NDim+i);
 
    for (int i=0; i<NDim; i++) {
      dstshape(i) *= appGridScale;
      dstextent(i) *= appGridScale;
    } 
 
    float dk=max(dstextent/dstshape);
    float kwidth=kernelwidth*dk*appGridScale;
    ODINLOG(odinlog,normalDebug) << "dstshape/dstextent/dk/kwidth=" << dstshape << "/" << dstextent << "/" << dk << "/" << kwidth << STD_endl;
 
    Data<float,NDim> sampldens(dstshape);
    sampldens=gridder.init(dstshape, dstextent, srccoords, kernel, kwidth);
//    sampldens.autowrite("sampldens.jdx");
 
 
 
    // Create magnitude correction for regridded data according to the gridding kernel used
    int deapo_os_factor=4; // oversampling to account for small kernels
    TinyVector<int,NDim> shape_os(deapo_os_factor*dstshape);
    ComplexData<NDim> kspace_kernel_os(shape_os);
    kspace_kernel_os=STD_complex(0.0);
 
    // fill k-space with gridding kernel
    for(unsigned int i=0; i<kspace_kernel_os.size(); i++) {
      TinyVector<int,NDim> index=kspace_kernel_os.create_index(i);
      TinyVector<int,NDim> dist=index-shape_os/2;
      float radius=secureDivision( sqrt(double(sum(dist*dist))), 0.5*deapo_os_factor*kernelwidth*appGridScale);
      if(radius<=1.0) kspace_kernel_os(index)=STD_complex(kernel.calculate(radius));
    }
 
 
    ODINLOG(odinlog,normalDebug) << "deapo_shift=" << deapo_shift << STD_endl;
 
    if(deapo_shift) {
      // Take spatial offset into account
      TinyVector<float,3> reloffset=controller.reloffset();
      TinyVector<float,NDim> reloffset_os;
      for(int i=0; i<NDim; i++) reloffset_os(i)=reloffset(3-NDim+i)/deapo_os_factor/appGridScale;
      if(max(abs(reloffset))>0.0) kspace_kernel_os.modulate_offset(reloffset_os);
    }
 
    // FFT and cutting out central part with de-apodize function
    kspace_kernel_os.fft();
    Data<float,NDim> apodize(dstshape);
    apodize=cabs(kspace_kernel_os(RectDomain<NDim>((shape_os-dstshape)/2, (shape_os-dstshape)/2+dstshape-1)));
//    apodize.autowrite("apodize"+ptos(traj)+".jdx");
 
    if(min(apodize)>0.0) {
      RecoData rdposted;
 
      // post separately for templates
      rdposted.coord().index[templtype]=coord.index[templtype];
      rdposted.coord().index[templtype].set_mode(RecoIndex::separate);
 
      rdposted.data(Rank<NDim>()).reference(float2real(mean(apodize)/apodize));
      controller.post_data(deapodize_post_str+postfix(), rdposted);
    } else {
      ODINLOG(odinlog,warningLog) << "FT of gridding kernel contains zero, not applying de-apodize correction" << STD_endl;
    }
 
 
    return gridder;
  }
  return it->second;
}
 
 
 
template<int NDim>
bool RecoGrid<NDim>::process(RecoData& rd, RecoController& controller) {
  Log<Reco> odinlog(c_label(),"process");
 
  if(!rd.coord().kspaceTraj) {
    ODINLOG(odinlog,errorLog) << "No kspaceTraj available for coord=" << rd.coord().print() << STD_endl;
    return false;
  }
 
  rd.coord().gridScaleFactor = gridScale;
 
  const Gridding<STD_complex,NDim>& gridder=get_gridder(rd.coord(), controller);
 
  TinyVector<int,NDim> outshape=gridshape(rd.coord(), controller);
  if (NDim > 1) {
    for (int i=0; i<NDim; i++) outshape(i) *= gridScale;
  }
 
  ComplexData<NDim> outdata(outshape);
 
 
  ComplexData<1>& adc=rd.data(Rank<1>());
  int adcSize=adc.extent(0);
  ODINLOG(odinlog,normalDebug) << "adcSize=" << adcSize << STD_endl;
 
  unsigned int offset=0;
  if(rd.coord().kspaceTraj->extent(0)>1) offset=rd.coord().index[cycle]*adcSize; // only for segmented 2D spirals, ignore cycle otherwise
 
  outdata=gridder(adc, offset);
 
  rd.data(Rank<NDim>()).reference(outdata);
 
  return execute_next_step(rd,controller);
}
 
 
template<int NDim>
bool RecoGrid<NDim>::query(RecoQueryContext& context) {
  if(context.mode==RecoQueryContext::prep) {
    context.controller.announce_data(deapodize_post_str+postfix());
  }
  return RecoStep::query(context);
}
 
 
 
template<int NDim>
bool RecoDeapodize<NDim>::process(RecoData& rd, RecoController& controller) {
 
  RecoData rddeapo;
  if(!controller.inquire_data(*this, deapodize_post_str+postfix(), rddeapo)) return false;
 
  ComplexData<3>& data=rd.data(Rank<3>());
  ComplexData<NDim>& deapo=rddeapo.data(Rank<NDim>());
 
  for(unsigned int i=0; i<data.size(); i++) {
    TinyVector<int,3> index=data.create_index(i);
 
    TinyVector<int,NDim> deapoindex;
    for(int i=0; i<NDim; i++) deapoindex(i)=index(3-NDim+i);
 
    data(index)*=deapo(deapoindex);
 
  }
 
 
//  for(int iphase3d=0; iphase3d<data.extent(0); iphase3d++) data(iphase3d,all,all)=data(iphase3d,all,all)*deapo(all,all);
 
  return execute_next_step(rd,controller);
}
 
 
 
 
bool RecoGridCut::process(RecoData& rd, RecoController& controller) {
 
  float scaleFactor = rd.coord().gridScaleFactor;
  if (scaleFactor != 1.0) {
 
    Range all=Range::all();
 
    ComplexData<3>& indata=rd.data(Rank<3>());
    TinyVector<int,3> inshape(indata.shape());
    TinyVector<int,3> outshape(inshape(0), inshape(1) / scaleFactor, inshape(2) / scaleFactor);
 
    TinyVector<int,2> phasebounds((inshape(1) - outshape(1)) / 2, (inshape(1) + outshape(1)) / 2);
    TinyVector<int,2> readbounds((inshape(1) - outshape(1)) / 2, (inshape(1) + outshape(1)) / 2);
 
    ComplexData<3> outdata(outshape);
    outdata = indata(all,Range(phasebounds(0),phasebounds(1)),Range(readbounds(0),readbounds(1)));
    rd.data(Rank<3>()).reference(outdata);
 
    rd.coord().gridScaleFactor = 1.0; // reset to 1.0 since the scaling has been removed from the data
  }
 
  return execute_next_step(rd,controller);
}