odin 1.8.4-1ubuntu2 / usr / share / odin / sequences / odinquant.cpp

#include <odinseq/seqall.h>


#define T2_CPMG_ECHO_FACTOR 2.5
#define T1_IR_FACTOR 1.2

class METHOD_CLASS : public SeqMethod {

 public:
  METHOD_CLASS(const STD_string& label);

  void method_pars_init();
  void method_seq_init();
  void method_rels();
  void method_pars_set();
  unsigned int numof_testcases() const {return 4;}

 private:

  JDXenum   Mode;
  JDXdouble ExpectedT1;
  JDXdouble ExpectedT2;
  JDXdouble MinEchoTime;
  JDXdouble MaxEchoTime;
  JDXint    NumOfEchoes;
  JDXint    NumOfDummyEchoes;
  JDXdoubleArr EchoTimes;
  JDXdouble B1PulseDur;

  SeqPulsar exc;
  SeqPulsar refoc;
  SeqPulsar inv;
  SeqAcqRead acqread;
  SeqDelay acqdummy;
  SeqAcqDeph readdeph;
  SeqObjLoop echoloop;
  SeqObjLoop dummyecholoop;
  SeqObjLoop peloop;
  SeqGradConstPulse spoiler1;
  SeqGradVectorPulse spoiler2;

  // objects for fieldmap scan
  SeqAcqEPI epi;
  SeqAcqDeph epideph;

  SeqVector echovec; // index vector for echoes

  SeqDelay excdelay;
  SeqDelay echodelay1;
  SeqDelay echodelay2;
  SeqDelay invdelay;
  SeqGradPhaseEnc pe1;
  SeqGradPhaseEnc pe2;
  SeqObjList echopart;
  SeqObjList grechpart;
  SeqObjList dummyechopart;
  SeqObjList kernel;
  SeqObjList slicepart;
  SeqDelay relaxdelay;
  SeqDelay dummype;
  SeqObjLoop sliceloop;
  SeqObjLoop acculoop;
  SeqObjLoop fliploop;

  SeqObjList postexc;
  SeqObjList postrefoc;
  SeqObjList postrefoctmpl;
  SeqObjList postinv;
  SeqObjList baselinepart;


  SeqGradEcho grech;

  fvector b1flips;

};

//////////////////////////////////////////////////////////////////////////////////////////

METHOD_CLASS::METHOD_CLASS (const STD_string& label)
                                   : SeqMethod(label) {
  set_description("A method to measure proton density, T1, T2, and frequency offset. "
                  "T1 maps and proton density maps are acquired using a Look-Locher sequence. "
                  "T2 maps are acquired using a CPMG sequence whereby each echo train acquires one line in k-space for all TEs. "
                  "Field maps are acquired using a spoiled gradient-echo sequence with different TEs. "
                  "B1 maps are acquired as described in Dowell et al., Magn Reson Med 58 (2007). "
                  "Maps are generated automatically by the reco method in this sequence module. ");
}


void METHOD_CLASS::method_pars_init() {
  Log<Seq> odinlog(this,"method_pars_init");

  Mode.add_item("T2_CPMG");
  Mode.add_item("T1_LookLocker");
  Mode.add_item("FIELD_MAP");
  Mode.add_item("B1");
  Mode.set_actual(get_current_testcase()); // alternative default settings for sequence test
  Mode.set_description("The quantitiy (and sequence) to measure.");

  ExpectedT1.set_minmaxval(10.0,10000.0);
  ExpectedT1=1300.0;
  ExpectedT1.set_description("The T1 which is expected. The TR and inversion-recovery timing will be optimized for this value.");
  ExpectedT1.set_unit(ODIN_TIME_UNIT);

  ExpectedT2.set_minmaxval(10.0,10000.0);
  ExpectedT2=90.0;
  ExpectedT2.set_description("The T2 which is expected. The CPMG timing will be optimized for this value.");
  ExpectedT2.set_unit(ODIN_TIME_UNIT);

  commonPars->set_MatrixSize(readDirection,128);
  commonPars->set_FlipAngle(90.0);


  NumOfEchoes=8;
  NumOfDummyEchoes=1;
  
  B1PulseDur=6.4; // use 30.0 for 7 Tesla;

  EchoTimes.resize(NumOfEchoes);

  append_parameter(Mode,"Mode");

  append_parameter(ExpectedT1,"ExpectedT1");
  append_parameter(ExpectedT2,"ExpectedT2");

  append_parameter(MinEchoTime,"MinEchoTime",noedit);
  append_parameter(MaxEchoTime,"MaxEchoTime",noedit);
  append_parameter(NumOfEchoes,"NumOfEchoes");
  append_parameter(NumOfDummyEchoes,"NumOfDummyEchoes");

  append_parameter(EchoTimes,"EchoTimes");

  append_parameter(B1PulseDur,"B1PulseDur");


}



void METHOD_CLASS::method_seq_init() {
  Log<Seq> odinlog(this,"method_seq_init");

  if(NumOfEchoes<2) NumOfEchoes=2;

  // adjust echo timings
  if(Mode=="T2_CPMG") {
    double te=ExpectedT2*T2_CPMG_ECHO_FACTOR/NumOfEchoes;
    commonPars->set_EchoTime(te,noedit);
    MinEchoTime=te;
    MaxEchoTime=NumOfEchoes*te;
  }
  if(Mode=="T1_LookLocker") {
    double ti_step=ExpectedT1*T1_IR_FACTOR/(NumOfEchoes-1);
    commonPars->set_EchoTime(ti_step,noedit);
    MinEchoTime=ExpectedT1/50.0;
    MaxEchoTime=MinEchoTime+(NumOfEchoes-1)*ti_step;
    float trwait=4.0*ExpectedT1;
    if(commonPars->get_RepetitionTime()<trwait) commonPars->set_RepetitionTime(trwait);
  }
  if(Mode=="FIELD_MAP") {
    if(NumOfEchoes%2) NumOfEchoes++; // even number of echoes (echo pairs) in EPI module
    // get TE later from EPI module
  }

  ///////////////// Pulses: /////////////////////

  float slicethick=geometryInfo->get_sliceThickness();
  float spatres=slicethick/4.0;

  float gamma=systemInfo->get_gamma();

  // Sinc Excitation Pulse
  exc=SeqPulsarSinc("exc",slicethick,true,4.0,commonPars->get_FlipAngle(),spatres,256);
  exc.set_freqlist( gamma * exc.get_strength() / (2.0*PII) *  geometryInfo->get_sliceOffsetVector() );
  exc.set_pulse_type(excitation);
  if(Mode=="T1_LookLocker") exc.set_rephased(false);

  if(Mode=="B1") {
    exc.set_rephased(false);
    b1flips.resize(5);
    b1flips[0]=110.0;  
    b1flips[1]=145.0;
    b1flips[2]=180.0;
    b1flips[3]=215.0;
    b1flips[4]=250.0;

    exc.set_flipangles(b1flips);
    exc.set_pulsduration(B1PulseDur);
  }

  // 90-180-90 Refocusing Pulse
  refoc=SeqPulsarBP("refoc", 1.0, 180.0);
  refoc.set_composite_pulse("90(x) 180(y) 90(x)");
  refoc.set_pulse_type(refocusing);


  // Inversion Pulse
  inv=SeqPulsar("inv",false,false);
  inv.set_dim_mode(zeroDeeMode);
  inv.set_Tp(10.0);
  inv.resize(256);
  inv.set_spat_resolution(0.0);
  inv.set_shape("Wurst(2.96,26.48)");  // bandwidth ~ 1kHz
  inv.set_trajectory("Const(0.0,1.0)");
  inv.set_filter("Gauss");
  inv.set_pulse_type(inversion);
  inv.refresh();




  ////////////////// Resolution: /////////////////////////////////

  // uniform resolution in read and phase direction
  float Resolution=secureDivision(geometryInfo->get_FOV(readDirection),commonPars->get_MatrixSize(readDirection));
  commonPars->set_MatrixSize(phaseDirection,int(secureDivision(geometryInfo->get_FOV(phaseDirection),Resolution)+0.5),noedit);


  //////////////// Phase Encoding: //////////////////////////

  // Phase encoding before acquisition
  pe1=SeqGradPhaseEnc("pe1",commonPars->get_MatrixSize(phaseDirection),geometryInfo->get_FOV(phaseDirection),
                      phaseDirection,0.25*systemInfo->get_max_grad(),
                      linearEncoding, noReorder, 1, commonPars->get_ReductionFactor(), DEFAULT_ACL_BANDS, commonPars->get_PartialFourier());



  // Rephasing before the refocusing pulse for CPMG
  pe2=pe1;
  pe2.set_label("pe2");
  pe2.invert_strength();


  SeqObjList pe1dummy; pe1dummy+=pe1;
  dummype=SeqDelay("dummype",pe1dummy.get_duration());


  //////////////// Readout: //////////////////////////////

  float os_read=1.0;

  acqread=SeqAcqRead("acqread",commonPars->get_AcqSweepWidth(),commonPars->get_MatrixSize(readDirection), geometryInfo->get_FOV(readDirection), readDirection, os_read);

  dephaseMode dephmode=spinEcho;
  if(Mode=="T1_LookLocker" || Mode=="FIELD_MAP") dephmode=FID;

  readdeph=SeqAcqDeph("readdeph",acqread,dephmode);

  acqdummy=SeqDelay("acqdummy",acqread.get_duration());


  //////////////// Gradient Echo Module: //////////////////////////////

  grech=SeqGradEcho("grech", exc, commonPars->get_AcqSweepWidth(),
              commonPars->get_MatrixSize(readDirection), geometryInfo->get_FOV(readDirection),
              commonPars->get_MatrixSize(phaseDirection), geometryInfo->get_FOV(phaseDirection),
              linearEncoding, noReorder, 1, commonPars->get_ReductionFactor(), DEFAULT_ACL_BANDS, false, commonPars->get_PartialFourier());

  //////////////// EPI module: //////////////////////////////

  epi=SeqAcqEPI("epi",commonPars->get_AcqSweepWidth(),
                commonPars->get_MatrixSize(readDirection), geometryInfo->get_FOV(readDirection),
                commonPars->get_MatrixSize(phaseDirection), geometryInfo->get_FOV(phaseDirection),
                commonPars->get_MatrixSize(phaseDirection), 1, os_read,"",0,0,linear,false,1.0,0.0,NumOfEchoes/2);

  epideph=SeqAcqDeph("epideph",epi,FID);


  //////////////// Timng Delays: //////////////////////////////

  excdelay=SeqDelay("excdelay",systemInfo->get_min_duration(delayObj));

  echodelay1=SeqDelay("echodelay1",systemInfo->get_min_duration(delayObj));

  echodelay2=SeqDelay("echodelay2",systemInfo->get_min_duration(delayObj));

  relaxdelay=SeqDelay("relaxdelay",0.0);

  invdelay=SeqDelay("invdelay");

  //////////////// Echo Vector: //////////////////////////////

  // an index vector is used which is responsible for assigning
  // the ADCs to the correct k-space data in the reconstruction

  echovec=SeqVector("echovec",NumOfEchoes);



  //////////////// Loops: //////////////////////////////

  echoloop=SeqObjLoop("echoloop");

  dummyecholoop=SeqObjLoop("dummyecholoop");

  peloop=SeqObjLoop("peloop");

  sliceloop=SeqObjLoop("sliceloop");

  acculoop=SeqObjLoop("acculoop");

  fliploop=SeqObjLoop("fliploop");


  //////////////// Spoiler: //////////////////////////////

  double SpoilerDuration=2.0;
  double SpoilerStrength=60.0;

  float spoiler_strength_phys=(double)SpoilerStrength/100.0*systemInfo->get_max_grad();
  spoiler1=SeqGradConstPulse("spoiler1",phaseDirection,spoiler_strength_phys,SpoilerDuration);


  fvector spoilertrims(NumOfEchoes);
  spoilertrims.fill_linear(-1.0,1.0);
  spoiler2=SeqGradVectorPulse("spoiler2",readDirection,spoiler_strength_phys,spoilertrims,SpoilerDuration*double(NumOfEchoes));



  //////////////// Building the sequence: //////////////////////////////

  postexc=SeqObjList("postexc");
  postrefoc=SeqObjList("postrefoc");
  postrefoctmpl=SeqObjList("postrefoctmpl");

  echopart=SeqObjList("echopart");
  grechpart=SeqObjList("grechpart");
  dummyechopart=SeqObjList("dummyechopart");
  kernel=SeqObjList("kernel");
  slicepart=SeqObjList("slicepart");

  if(Mode=="T2_CPMG") {

    postexc=excdelay + readdeph ;
    postrefoc=    spoiler1  + echodelay1 + pe1;
    postrefoctmpl=spoiler1  + echodelay1 + dummype;

    echopart=      spoiler1 + refoc + postrefoc     + acqread  + pe2     + echodelay2 ;
    dummyechopart= spoiler1 + refoc + postrefoctmpl + acqdummy + dummype + echodelay2 ;

    if (NumOfDummyEchoes>1) {
      kernel= exc + postexc + dummyecholoop(dummyechopart)[NumOfDummyEchoes] + echoloop(echopart)[echovec];
    } else {
      kernel= exc + postexc + echoloop(echopart)[echovec];
    }

    slicepart= sliceloop( kernel + relaxdelay )[exc];

    set_sequence( acculoop( peloop( slicepart )[pe1][pe2] )[commonPars->get_NumOfRepetitions()] );
  }


  if(Mode=="T1_LookLocker") {

   grechpart= sliceloop( grech )[exc];

   echopart=  grechpart + spoiler2 + echodelay1;

   kernel= inv + spoiler1 + invdelay + echoloop( echopart )[spoiler2];

   slicepart=  kernel + relaxdelay;

   set_sequence(  acculoop (  peloop( slicepart )[grech.get_pe_vector()] )[commonPars->get_NumOfRepetitions()] );
  }



  if(Mode=="FIELD_MAP") {

    kernel= exc + epideph + epi;

    slicepart= sliceloop(  kernel + relaxdelay  )[exc];

    set_sequence(  acculoop ( peloop( slicepart )[epideph.get_epi_segment_vector()] )[commonPars->get_NumOfRepetitions()] );
  }

  if(Mode=="B1") {

   slicepart= sliceloop(  grech + spoiler1 + relaxdelay  )[exc];

   set_sequence(
     fliploop(
       acculoop (
         peloop(
           slicepart
         )[grech.get_pe_vector()]
       )[commonPars->get_NumOfRepetitions()]
     )[exc.get_flipangle_vector()]
   );

  }


}


void METHOD_CLASS::method_rels() {
  Log<Seq> odinlog(this,"method_rels");

  //////////////// Timings: //////////////////////////////

  float exc2refoc_duration;
  float refoc2acq_duration;
  float acq_shift;

  float min_echo_time=0.0;

  float mindelaydur=systemInfo->get_min_duration(delayObj);

  echodelay1.set_duration(mindelaydur);
  echodelay2.set_duration(mindelaydur);
  excdelay.set_duration(mindelaydur);

  EchoTimes.resize(NumOfEchoes);



  if(Mode=="T2_CPMG") {

    //////////////// Duration from the middle of the excitation pulse ///////////////////
    //////////////// to the middle of the refocusing pulse: /////////////////////////////
    exc2refoc_duration=( exc.get_duration() - exc.get_magnetic_center() )
                          + postexc.get_duration()
                          + echoloop.get_preduration()
                          + spoiler1.get_duration()
                          + refoc.get_magnetic_center();



    if(min_echo_time < exc2refoc_duration*2.0 ) min_echo_time=exc2refoc_duration*2.0;

    float echoloopdur=echoloop.get_preduration()+echopart.get_duration();

    if(min_echo_time < echoloopdur) min_echo_time=echoloopdur;

    if(commonPars->get_EchoTime()<min_echo_time) commonPars->set_EchoTime(min_echo_time);

    echodelay1.set_duration( ( commonPars->get_EchoTime() - echoloopdur )/2.0 );
    echodelay2.set_duration( ( commonPars->get_EchoTime() - echoloopdur )/2.0 );

    //////////////// Duration from the middle of the refocusing pulse ///////////////////
    //////////////// to the middle of the acquisition window: ///////////////////////////
    refoc2acq_duration=( refoc.get_duration() - refoc.get_magnetic_center() )
                     + postrefoc.get_duration()
                     + acqread.get_acquisition_center();


    echoloopdur=echoloop.get_preduration()+echopart.get_duration();

    acq_shift=refoc2acq_duration-echoloopdur/2.0;
    echodelay1.set_duration( echodelay1.get_duration() - acq_shift );
    echodelay2.set_duration( echodelay2.get_duration() + acq_shift );


    excdelay=commonPars->get_EchoTime()/2.0-exc2refoc_duration;

    EchoTimes.fill_linear(commonPars->get_EchoTime(),NumOfEchoes*commonPars->get_EchoTime());
  }






  if(Mode=="T1_LookLocker") {

    float inv2acq_duration=( inv.get_duration() - inv.get_magnetic_center() )
                           +spoiler1.get_duration()
                           +grech.get_acquisition_center();

    if(MinEchoTime<inv2acq_duration) MinEchoTime=inv2acq_duration;

    float oneloopdur=echopart.get_duration();

    float lastechotime=MinEchoTime+(NumOfEchoes-1)*oneloopdur;

    if(MaxEchoTime<lastechotime) MaxEchoTime=lastechotime;

    float echotimediff=MaxEchoTime-MinEchoTime;

    echodelay1=echotimediff/float(NumOfEchoes-1)-oneloopdur;

    invdelay=MinEchoTime-inv2acq_duration;

    EchoTimes.fill_linear(MinEchoTime,MaxEchoTime);
  }

  if(Mode=="FIELD_MAP" || Mode=="B1") {
    MinEchoTime=0.0;
    MaxEchoTime=0.0;
    EchoTimes.fill_linear(0.0,0.0);
  }
  ODINLOG(odinlog,normalDebug) << "EchoTimes=" << EchoTimes.printbody() << STD_endl;


  relaxdelay=0.0;
  double slicepackdur=slicepart.get_duration();
  ODINLOG(odinlog,normalDebug) << "slicepackdur=" << slicepackdur << STD_endl;

  int n_slices_to_consider=geometryInfo->get_nSlices();
  if(Mode=="T1_LookLocker") n_slices_to_consider=1; // because of inversion pulse
  if(Mode=="T2_CPMG") n_slices_to_consider=1; // because of hard refoc pulse

  if(commonPars->get_RepetitionTime()<slicepackdur) commonPars->set_RepetitionTime(slicepackdur);
  else relaxdelay=secureDivision(commonPars->get_RepetitionTime()-slicepackdur,n_slices_to_consider);



  float flipangle=90.0;

  if(Mode=="FIELD_MAP") {
    // calculate Ernst angle according to TR
    flipangle=180.0/PII * acos( exp ( -secureDivision ( commonPars->get_RepetitionTime(), ExpectedT1) ) );
    commonPars->set_FlipAngle( flipangle, noedit );
  } else if(Mode=="T1_LookLocker") {
    flipangle=commonPars->get_FlipAngle();
    float maxflip=15.0;
    if(flipangle>maxflip) flipangle=maxflip;
    commonPars->set_FlipAngle( flipangle, edit );
  } else {
    commonPars->set_FlipAngle( flipangle, noedit );
  }

  if(Mode!="B1") {
    exc.set_flipangle( flipangle );
  }

}


void METHOD_CLASS::method_pars_set() {
  Log<Seq> odinlog(this,"method_pars_set");

  if(Mode=="T2_CPMG") {

    // inform the readout about the used phase encoding and slice vector (for automatic reconstruction)
    acqread.set_reco_vector(line,pe1).set_reco_vector(slice,exc);

    // store different echoes in te dimension with echo times attached
    acqread.set_reco_vector(userdef, echovec, EchoTimes);

    recoInfo->set_PostProc3D("usercoll | expfit(true) ");
  }


  if(Mode=="T1_LookLocker") {

    // store different echoes in te dimension with echo times attached
    grech.set_reco_vector(userdef, spoiler2, EchoTimes);

    recoInfo->set_PostProc3D("usercoll | t1fit ");
  }


  if(Mode=="FIELD_MAP") {

    // inform the readout about the used slice vector (for automatic reconstruction)
    epi.set_reco_vector(slice,exc);

//    epi.set_te_offset(exc.get_duration()-exc.get_magnetic_center()+epideph.get_duration());

    recoInfo->set_PreProc3D("tecoll | fieldmap "); // calculate fieldmap before summing up channels!
  }

  if(Mode=="B1") {

    // store different echoes in te dimension with echo times attached
    grech.set_reco_vector(userdef, exc.get_flipangle_vector(), fvector2dvector(b1flips));

    recoInfo->set_PostProc3D("usercoll | b1fit ");
  }

}


///////////////////////////////////////////////////////////////////////////////////////////


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT