grappa_code.h

#include "grappa.h"
#include "data.h"
#include "controller.h"
 
#include <odindata/linalg.h>
 
 
STD_string grappa_postlabel(recoDim dim) {
  return STD_string("grappaweights_")+recoDimLabel[dim];
}
 
 
 
ivector next_neighbours_offset(unsigned int reductionFactor, int numof_neighb, int ired) {
  Log<Reco> odinlog("","next_neighbours_offset");
 
  ivector result(numof_neighb);
 
  int negfirst=(ired<int(reductionFactor/2)); // determine initial direction
 
  for(int i=0; i<numof_neighb; i++) {
    int i2=i/2;  // interleaved positive and negative direction
    if((i+negfirst)%2) result[i]=-ired-1-i2*reductionFactor; // negative direction
    else               result[i]=reductionFactor-ired-1+i2*reductionFactor; // positive direction
  }
 
  ODINLOG(odinlog,normalDebug) << "result[" << reductionFactor << "/" << ired << "]" << result << STD_endl;
 
  return result;
}
 
 
bool measured_line(const ivector& indexvec, int index) {
  for(unsigned int i=0; i<indexvec.size(); i++) {
    if(indexvec[i]==index) return true;
  }
  return false;
}
 
 
void set_separate_index(RecoCoord& coord, recoDim orthoDim, int iortho) {
  coord.set_mode(RecoIndex::separate, orthoDim);
  coord.index[orthoDim]=iortho;
}
 
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex, class Ignore>
void RecoGrappaWeights<interpolDim,orthoDim,interpolIndex,orthoIndex,Ignore>::init() {
 
  reduction_factor=0;
  reduction_factor.set_description("Reduction factor");
  append_arg(reduction_factor,"reduction_factor");
 
  neighbours_read=5;
  neighbours_read.set_cmdline_option("gr").set_description("Number of neigbours in read direction used for GRAPPA interpolation");
  append_arg(neighbours_read,"neighbours_read");
 
  neighbours_phase=2;
  neighbours_phase.set_cmdline_option("gp").set_description("Number of neigbouring measured k-space lines used for GRAPPA interpolation");
  append_arg(neighbours_phase,"neighbours_phase");
 
  svd_trunc=0.001;
  svd_trunc.set_cmdline_option("gs").set_description("Truncation value of SVD (i.e. regularization) when calculating GRAPPA weights");
  append_arg(svd_trunc,"svd_trunc");
 
  discard_level=0.0; //0.01;
  discard_level.set_cmdline_option("gd").set_description("Fraction of k-psace points to discard in auto-calibration lines because of high residuals");
  append_arg(discard_level,"discard_level");
}
 
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex, class Ignore>
bool RecoGrappaWeights<interpolDim,orthoDim,interpolIndex,orthoIndex,Ignore>::process(RecoData& rd, RecoController& controller) {
  Log<Reco> odinlog(c_label(),"process");
 
  RecoCoord coordcopy(rd.coord()); // local copy to modify coord according to ignored dims
  Ignore::modify(RecoIndex::ignore, coordcopy);
 
  RecoData rdweights(coordcopy);
  if(!calc_weights(rdweights.data(Rank<5>()), coordcopy, rd.data(Rank<4>()))) return false;
 
  controller.post_data(grappa_postlabel(interpolDim), rdweights);
  ODINLOG(odinlog,normalDebug) << "Posted rdweights=" << rdweights.coord().print() << STD_endl;
 
  return execute_next_step(rd,controller);
}
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex, class Ignore>
bool RecoGrappaWeights<interpolDim,orthoDim,interpolIndex,orthoIndex,Ignore>::query(RecoQueryContext& context) {
  Log<Reco> odinlog(c_label(),"query");
  if(context.mode==RecoQueryContext::prep) {
    context.controller.announce_data(grappa_postlabel(interpolDim));
    RecoCoord coordcopy(context.coord);
    Ignore::modify(RecoIndex::ignore, coordcopy);
    if(!measlines.init(coordcopy, context.controller)) return false;
  }
  return RecoStep::query(context);
}
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex, class Ignore>
bool RecoGrappaWeights<interpolDim,orthoDim,interpolIndex,orthoIndex,Ignore>::calc_weights(ComplexData<5>& weights, const RecoCoord& trainingcoord, const ComplexData<4>& trainingdata) {
  Log<Reco> odinlog(c_label(),"calc_weights");
 
  Range all=Range::all();
 
  ODINLOG(odinlog,normalDebug) << "neighbours_read/neighbours_phase=" << neighbours_read << "/" << neighbours_phase << STD_endl;
 
  TinyVector<int,4> inshape=trainingdata.shape();
  int nChannels=inshape(0);
  int sizeOrtho=inshape(1+orthoIndex);
  int sizeInterpol=inshape(1+interpolIndex);
  int sizeRead=inshape(3);
  ODINLOG(odinlog,normalDebug) << "nChannels/sizeOrtho/sizeInterpol/sizeRead/reduction_factor=" << nChannels << "/" << sizeOrtho << "/" << sizeInterpol << "/" << sizeRead << "/" << reduction_factor << STD_endl;
 
  // Shape for 3D volumes
  TinyVector<int,3> shape3d(inshape(1),inshape(2),inshape(3));
 
  // shape to calculate index vector within interpolation net
  TinyVector<int,3> netshape(nChannels, neighbours_phase, neighbours_read);
 
  int ncols=product(netshape); // Size of the interpolation net
  ComplexData<1> weights_net_vec(ncols); // Interpolation net around a particular root node
 
  weights.resize(nChannels, reduction_factor-1, nChannels, neighbours_phase, neighbours_read); // the result
  weights=STD_complex(0.0);
 
  for(int ired=0; ired<(reduction_factor-1); ired++) { // loop over (R-1) missing lines 
 
    // Mask for the root nodes available to calculate the weights
    Data<int,4> acl_mask(inshape);
    acl_mask=0;
 
    RecoCoord aclcoord(trainingcoord);
    for(int iortho=0; iortho<sizeOrtho; iortho++) {
 
       // GRAPPA encoding might be different in each of the orthogonal phase encoding partitions
      set_separate_index(aclcoord, orthoDim, iortho);
 
      ivector aclindices=acl_lines(aclcoord, neighbours_phase, ired);
      ODINLOG(odinlog,normalDebug) << "aclindices(" << ired << ", " << aclcoord.print() << ")=" << aclindices << STD_endl;
 
      for(int iacl=0; iacl<int(aclindices.size()); iacl++) {
        TinyVector<int,2> aclindex;
        aclindex(orthoIndex)=iortho;
        aclindex(interpolIndex)= aclindices[iacl];
        acl_mask(all, aclindex(0), aclindex(1), Range( neighbours_read/2, sizeRead-1-neighbours_read/2 ))=1;
      }
    }
 
    TinyVector<int,3> acl_signal_shape(shape3d); // used only to calculate index vector of root nodes
 
    // Offset of neighbouring lines used to calculate weights for this particular reduction index
    ivector neighboffset=next_neighbours_offset(reduction_factor, neighbours_phase, ired);
    ODINLOG(odinlog,normalDebug) << "neighboffset(" << ired << ")=" << neighboffset.printbody() << STD_endl;
 
    for(int ichan_dst=0; ichan_dst<nChannels; ichan_dst++) { // calculate weights separately for each channel
 
      int nrows=sum(acl_mask(ichan_dst,all,all,all));
      ODINLOG(odinlog,normalDebug) << "ncols/nrows=" << ncols << "/" << nrows << STD_endl;
      if(ncols>nrows) {
        ODINLOG(odinlog,errorLog) << "Not enough auto-calibration data for " << netshape << " interpolation net" << STD_endl;
        return false;
      }
 
      // fill acl_signal_vec with signal values from ACLs
      ComplexData<1> acl_signal_vec(nrows); // Will hold root node data of a particular channel
      Array<TinyVector<int,3>,1> maskindex(nrows); // cache for coordinates of signal values contributing to SVD
      int irow=0;
      for(int i=0; i<product(acl_signal_shape); i++) {
        TinyVector<int,3> rowindex=index2extent(acl_signal_shape, i);
        TinyVector<int,4> aclindex(ichan_dst, rowindex(0), rowindex(1), rowindex(2));
        if(acl_mask(aclindex)) {
          maskindex(irow)=rowindex;
          STD_complex aclval=trainingdata(aclindex);
          if(cabs(aclval)==0.0) {
            ODINLOG(odinlog,warningLog) << "Zero acl at " << aclindex << STD_endl;
          }
          acl_signal_vec(irow)=aclval;
          irow++;
        }
      }
 
      // fill Matrix for SVD
      ComplexData<2> Matrix(nrows,ncols);
      for(int irow=0; irow<nrows; irow++) { // loop over root nodes
        TinyVector<int,3> aclindex=maskindex(irow);
 
        for(int icol=0; icol<ncols; icol++) {  // loop over interpolation net, i.e. over neighbourhood of root node
          TinyVector<int,3> windex=index2extent(netshape,icol);
 
          int ichan_src=windex(0);
          int ipoloffset=neighboffset[windex(1)];
          int readoffset=windex(2)-neighbours_read/2; // symmetrical about root node
 
          TinyVector<int,4> trainingindex;
          trainingindex(0)=ichan_src;
          trainingindex(1)=aclindex(0);
          trainingindex(2)=aclindex(1);
          trainingindex(3)=aclindex(2)+readoffset;
 
          trainingindex(1+interpolIndex)+=ipoloffset; // take neighbour in interpolation direction
 
          STD_complex trainingval=trainingdata(trainingindex);
          if(cabs(trainingval)==0.0) {
            ODINLOG(odinlog,warningLog) << "Zero trainingdata at " << trainingindex << STD_endl;
          }
 
          Matrix(irow,icol)=trainingval;
        }
      }
 
      // solve system of linear equations using complex SVD
      weights_net_vec=solve_linear(Matrix, acl_signal_vec, svd_trunc);
 
 
      // Use residuals to eliminate noisy ACLs (Huo et al., JMRI 2008, 27:1412)
      if(discard_level>0.0) {
        Data<float,1> residuals(cabs(matrix_product(Matrix, weights_net_vec)-acl_signal_vec));
        int ndiscard=int(discard_level*nrows+0.5);
        if(ncols>(nrows-ndiscard)) {
          ODINLOG(odinlog,warningLog) << "Limiting ndiscard to " << ndiscard << " for sufficient auto-calibration data" << STD_endl;
          ndiscard=nrows-ncols;
        }
        for(int i=0; i<ndiscard; i++) {
          int irow_max=maxIndex(residuals)(0);
          acl_signal_vec(irow_max)=STD_complex(0.0);
          Matrix(irow_max,all)=STD_complex(0.0);
          residuals(irow_max)=0.0;
        }
 
        // solve again
        weights_net_vec=solve_linear(Matrix, acl_signal_vec, svd_trunc);
      }
 
      for(int icol=0; icol<ncols; icol++) {
        TinyVector<int,3> windex=index2extent(netshape,icol);
        weights(ichan_dst, ired, windex(0), windex(1), windex(2))=weights_net_vec(icol);
      }
 
    } // end loop over channels
 
  } // end loop over reduction inidices
 
  ODINLOG(odinlog,normalDebug) << "cabs(sum(weights" << trainingcoord.print() << "))=" << cabs(sum(weights)) << STD_endl;
 
  return true;
}
 
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex, class Ignore>
ivector RecoGrappaWeights<interpolDim,orthoDim,interpolIndex,orthoIndex,Ignore>::acl_lines(const RecoCoord& coord, int numof_neighb, int ired) const {
  Log<Reco> odinlog(c_label(),"acl_lines");
 
  ivector neighboffset=next_neighbours_offset(reduction_factor, numof_neighb, ired);
  ODINLOG(odinlog,normalDebug) << "neighboffset=" << neighboffset.printbody() << STD_endl;
 
  ivector indexvec=measlines.get_indices(coord);
  ODINLOG(odinlog,normalDebug) << "indexvec(" << coord.print() << ")=" << indexvec.printbody() << STD_endl;
 
  STD_list<int> indexlist;
  for(unsigned int i=0; i<indexvec.size(); i++) {
    int iline=indexvec[i];
    bool acl=true;
    for(int ineighb=0; ineighb<numof_neighb; ineighb++) { // check whether all neighbours are there
      int neighbindex=iline+neighboffset[ineighb];
      if(!measured_line( indexvec, neighbindex)) acl=false;
    }
    if(acl) indexlist.push_back(iline);
  }
 
  return list2vector(indexlist);
}
 
 
 
 
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex>
void RecoGrappa<interpolDim,orthoDim,interpolIndex,orthoIndex>::init() {
 
  reduction_factor=0;
  reduction_factor.set_description("Reduction factor");
  append_arg(reduction_factor,"reduction_factor");
 
  keep_shape=false;
  keep_shape.set_description("Do not correct shape of interpolated k-space according to image space size");
  append_arg(keep_shape,"keep_shape");
 
}
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex>
bool RecoGrappa<interpolDim,orthoDim,interpolIndex,orthoIndex>::process(RecoData& rd, RecoController& controller) {
  Log<Reco> odinlog(c_label(),"process");
 
  if(!controller.data_announced(grappa_postlabel(interpolDim))) {
    ODINLOG(odinlog,errorLog) << "GRAPPA weights not available" << STD_endl;
    return false;
  }
 
  ODINLOG(odinlog,normalDebug) << "reduction_factor=" << reduction_factor << STD_endl;
 
  if(!keep_shape) {
    if(!correct_shape(rd, controller)) return false; // adjust shape to account for incomplete sampling at end of k-space
  }
 
  ODINLOG(odinlog,normalDebug) << "Requesting weights for " << rd.coord().print() << STD_endl;
  RecoData rdweights(rd.coord());
  if(!controller.inquire_data(*this, grappa_postlabel(interpolDim), rdweights)) return false;
  ODINLOG(odinlog,normalDebug) << "Retrieved weights with " << rdweights.coord().print() << STD_endl;
 
  const ComplexData<5>& weights=rdweights.data(Rank<5>());
  ODINLOG(odinlog,normalDebug) << "cabs(sum(weights(" << rdweights.coord().print() << ")))=" << cabs(sum(weights)) << STD_endl;
 
 
  ComplexData<4>& kspace=rd.data(Rank<4>());
  TinyVector<int,4> kspaceshape=kspace.shape();
  int nChannels=kspaceshape(0);
  int sizeOrtho=kspaceshape(1+orthoIndex);
  int sizeInterpol=kspaceshape(1+interpolIndex);
  int sizeRead=kspaceshape(3);
  ODINLOG(odinlog,normalDebug) << "nChannels/sizeOrtho/sizeInterpol/sizeRead/reduction_factor=" << nChannels << "/" << sizeOrtho << "/" << sizeInterpol << "/" << sizeRead << "/" << reduction_factor << STD_endl;
 
 
  int neighbours_phase=weights.extent(3);
  int neighbours_read=weights.extent(4);
  ODINLOG(odinlog,normalDebug) << "neighbours_phase/neighbours_read=" << neighbours_phase << "/" << neighbours_read << STD_endl;
 
  CoordCountMap interpolated_coords;
  if(rd.interpolated) {
    measlines.update(*(rd.interpolated)); // take previously interpolated coords into account
    interpolated_coords=(*(rd.interpolated)); // merge previously interpolated coords
  }
 
  // shape to calculate index vector within interpolation net
  TinyVector<int,3> netshape(nChannels, neighbours_phase, neighbours_read);
 
  // separate array for interpolated values to avoid recursive summation
  ComplexData<4> kspace_interpol(kspaceshape);
  kspace_interpol=STD_complex(0.0);
 
 
  for(int iortho=0; iortho<sizeOrtho; iortho++) {
 
    // GRAPPA encoding might be different in each of the orthogonal phase encoding partitions
    RecoCoord meascoord(rd.coord());
    set_separate_index(meascoord, orthoDim, iortho);
 
    ivector indexvec=measlines.get_indices(meascoord);
    ODINLOG(odinlog,normalDebug) << "indexvec(" << meascoord.print() << ")=" << indexvec.printbody() << STD_endl;
 
    if(indexvec.size()) {  // discard empty lines in orthogonal direction
 
      for(int ipol=0; ipol<sizeInterpol; ipol++) {
 
        if(!measured_line(indexvec, ipol)) {
 
          int ired=reduction_index(indexvec, sizeInterpol, ipol);
          ODINLOG(odinlog,normalDebug) << "ired(" << ipol << ")=" << ired << STD_endl;
 
          if(ired>=0 && ired<(reduction_factor-1)) { // omit lines which are discarded due to partial Fourier
 
            set_separate_index(meascoord, interpolDim, ipol);
            interpolated_coords[meascoord]++;
            ODINLOG(odinlog,normalDebug) << "Interpolating coord " << meascoord.print() << STD_endl;
 
            ivector neighboffset=next_neighbours_offset(reduction_factor, neighbours_phase, ired);
            ODINLOG(odinlog,normalDebug) << "neighboffset" << neighboffset << STD_endl;
 
            for(int iread=neighbours_read/2; iread<(sizeRead-neighbours_read/2); iread++) { // avoid interpolation of edges without neighbours
              for(int ichan=0; ichan<nChannels; ichan++) {
 
                TinyVector<int,4> kspaceindex;
                kspaceindex(0)=ichan;
                kspaceindex(1+orthoIndex)=iortho;
                kspaceindex(1+interpolIndex)=ipol;
                kspaceindex(3)=iread;
 
                STD_complex interpolval(0.0);
 
                // iterate over 3D interpolation net to accumulate signal value at root node
                for(int i=0; i<product(netshape); i++) {
                  TinyVector<int,3> netindex=index2extent(netshape,i); // index within net
                  int jchan=netindex(0);
                  int jphase=netindex(1);
                  int jread=netindex(2);
 
                  int offset_ipol=neighboffset[jphase];
                  int offset_read=jread-neighbours_read/2;  // symmetrical about root node
 
                  int src_ipol  =ipol+offset_ipol;
                  int src_iread =iread+offset_read;
 
                  if(src_ipol>=0 && src_ipol<sizeInterpol && src_iread>=0 && src_iread<sizeRead) { // check if src index is within k-space
 
                    TinyVector<int,4> srcindex;
                    srcindex(0)=jchan;
                    srcindex(1+orthoIndex)=iortho; // Use only src points from the same partition in orthogonal direction
                    srcindex(1+interpolIndex)=src_ipol;
                    srcindex(3)=src_iread;
 
                    STD_complex srcval=kspace(srcindex);
/*
                    if(cabs(srcval)==0.0) {
                      ODINLOG(odinlog,warningLog) << "Zero srcval at " << srcindex << " while interpolating " << kspaceindex << STD_endl;
                    }
*/
 
                    interpolval += srcval * weights(ichan,ired,jchan,jphase,jread);  // linear interpolation of kspace
                  }
                }
 
                kspace_interpol(kspaceindex)=interpolval;
              }
            }
          }
        }
      }
    }
  }
 
  kspace=kspace+kspace_interpol;
 
  rd.interpolated=&interpolated_coords; // Inform subsequent steps about the interpolated coordinates
 
  return execute_next_step(rd,controller);
}
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex>
bool RecoGrappa<interpolDim,orthoDim,interpolIndex,orthoIndex>::correct_shape(RecoData& rdkspace, const RecoController& controller) const {
  Log<Reco> odinlog(c_label(),"correct_shape");
 
  Range all=Range::all();
 
  ComplexData<4>& inkspace=rdkspace.data(Rank<4>());
  TinyVector<int,4> inshape=inkspace.shape();
 
  int nlines_dst=controller.image_size()(interpolIndex);
  int nlines_src=inshape(1+interpolIndex);
  if(nlines_dst==nlines_src) return true;
 
  if(nlines_dst<nlines_src) nlines_dst=nlines_src;
  TinyVector<int,4> outshape(inshape);
  outshape(1+interpolIndex)=nlines_dst;
 
  ODINLOG(odinlog,normalDebug) << "inshape/outshape=" << inshape << "/" << outshape << STD_endl;
 
  ComplexData<4> outkspace(outshape);
  outkspace=STD_complex(0.0);
 
  Range srcrange(0,nlines_src-1);
  if(interpolDim==line)   outkspace(all,all,srcrange,all)=inkspace(all,all,srcrange,all);
  if(interpolDim==line3d) outkspace(all,srcrange,all,all)=inkspace(all,srcrange,all,all);
  rdkspace.data(Rank<4>()).reference(outkspace);
 
  return true;
}
 
 
 
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex>
bool RecoGrappa<interpolDim,orthoDim,interpolIndex,orthoIndex>::query(RecoQueryContext& context) {
  Log<Reco> odinlog(c_label(),"query");
  if(context.mode==RecoQueryContext::prep) {
    if(!measlines.init(context.coord, context.controller)) return false;
  }
 
  return RecoStep::query(context);
}
 
 
template<recoDim interpolDim, recoDim orthoDim, int interpolIndex, int orthoIndex>
int RecoGrappa<interpolDim,orthoDim,interpolIndex,orthoIndex>::reduction_index(const ivector& indexvec, int sizePhase, int iphase) const {
  Log<Reco> odinlog(c_label(),"reduction_index");
 
  ODINLOG(odinlog,normalDebug) << "reduction_factor/sizePhase/iphase=" << reduction_factor << "/" << sizePhase << "/" << iphase << STD_endl;
 
  int posnext=-1;
  int negnext=-1;
 
  // search for next scanned line in positive direction
  for(int i=iphase; i<sizePhase; i++) {
    if(measured_line(indexvec, i)) {
      posnext=i-iphase;
      break;
    }
  }
 
  // search for next scanned line in negative direction
  for(int i=iphase; i>=0; i--) {
    if(measured_line(indexvec, i)) {
      negnext=iphase-i;
      break;
    }
  }
 
  ODINLOG(odinlog,normalDebug) << "posnext/negnext(" << iphase << ")=" << posnext << "/" << negnext << STD_endl;
 
  if(posnext<0 || negnext<0) return -1; // Measured line is missing in one or more directions -> Line is not surrounded by measured lines
 
  return negnext-1;
}