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

#include <odinseq/seqall.h>

class METHOD_CLASS : public SeqMethod {

private:
  JDXint    Blades;
  JDXbool   ShortAxis;
  JDXfloat  BladeOversampling;
  JDXbool   TakeMinEchoTime;
  JDXbool   RampSampling;
  JDXenum   RampMode;
  JDXfloat  RampSteepness;
  JDXbool   FatSaturation;
  JDXfloat  T1Ernst;
  JDXint    DummyCycles;
  JDXint    NumOfGradEchoes;
  JDXint    NumOfSamples;
  JDXdouble PulseDur;
  JDXbool   fMRITrigger;
  JDXbool   UpDownBlade;
  JDXbool   FieldMap;
  JDXbool   EPIRef;


  SeqPulsar  exc;
  SeqSat     fatsat;
  SeqAcqEPI  epiacq;
  SeqDelay   epiacq_dummy;
  SeqAcqEPI  epiacq_template;
  SeqAcqDeph deph;
  SeqAcqDeph deph_template;
  SeqObjLoop sliceloop;
  SeqObjLoop reploop;
  SeqObjLoop bladeloop;
  SeqObjLoop dummyloop;
  SeqDelay   trdelay;
  SeqObjList scan;
  SeqObjList dummypart;
  SeqObjList templatepart;
  SeqObjList imagingpart;
  SeqObjList slicepart;
  SeqObjList slicepart_dummy;
  SeqObjList slicepart_template;
  SeqObjList preppart;
  SeqDelay   exc2acq;
  SeqTrigger trigger;
  SeqGradTrapezParallel crusher;
  SeqDelay   crusherdelay;
  SeqFieldMap fmapscan;
  SeqVecIter phaseiter;

  SeqRotMatrixVector bladerot;
  SeqMagnReset reset;


  // stuff for EPIRef
  SeqAcqEPI  epiacq_ref;
  SeqDelay   epiacq_ref_dummy;
  SeqAcqEPI  epiacq_ref_template;
  SeqAcqEPI  epiacq_ref_grappa;
  SeqAcqDeph deph_ref;
  SeqAcqDeph deph_ref_template;
  SeqAcqDeph deph_ref_grappa;
  SeqObjLoop grappaloop_ref;
  SeqDelay   trdelay_ref;
  SeqObjList templatepart_ref;
  SeqObjList grappapart_ref;
  SeqObjList imagingpart_ref;
  SeqObjList slicepart_grappa_ref;
  SeqObjList slicepart_template_ref;
  SeqDelay   exc2acq_ref;


  dvector bladeangels;




public:

// This constructor creates an empty EPI sequence
METHOD_CLASS(const STD_string& label) : SeqMethod(label) {
  set_description("2D-PROPELLER-EPI with EPI reference scan.");
}


void method_pars_init() {

  // In this function, parameters are initialized and default values are set

  commonPars->set_MatrixSize(readDirection,128);
  commonPars->set_MatrixSize(phaseDirection,128,noedit);
  commonPars->set_NumOfRepetitions(1);
  commonPars->set_RepetitionTime(1000.0);
  commonPars->set_AcqSweepWidth(100.0);


  Blades=8;
  Blades.set_description("Number of PROPELLER blades");

  ShortAxis=false;
  ShortAxis.set_description("Readout direction along short axis of blades");

  BladeOversampling=1.0;
  BladeOversampling.set_minmaxval(0.0,100.0);
  BladeOversampling.set_unit("%").set_description("Oversampling (overlap) in short direction of blade");


  RampSampling=false;
  RampSampling.set_description("Perform sampling during gradient ramps");

  RampMode.add_item("linear",linear);
  RampMode.add_item("sinusoidal",sinusoidal);
  RampMode.add_item("half_sinusoidal",half_sinusoidal);
  RampMode.set_actual(linear);
  RampMode.set_description("The shape of the ramps of the read gradient");

  RampSteepness=1.0;
  RampSteepness.set_description("Relative steepness (slew rate) of the EPI readout ramps");


  FatSaturation=true;
  FatSaturation.set_description("Saturation of fat resonance prior to excitation");

  T1Ernst=1300.0;
  T1Ernst.set_minmaxval(0.0,5000.0).set_description("If non-zero, the flip angle will be set to the Ernst angle using this T1 for optimum SNR");

  TakeMinEchoTime=true;
  TakeMinEchoTime.set_description("Use minimum possible TE");

  DummyCycles=3;
  DummyCycles.set_description("Number of dummy shots before actual acquisition");

  fMRITrigger=true;
  fMRITrigger.set_description("External triggering");

  UpDownBlade=false;
  UpDownBlade.set_description("Sample each blade twice with opposite phase-encoding direction");

  FieldMap=false;
  FieldMap.set_description("Fieldmap pre-scan for distortion correction");

  PulseDur=2.0; // avoid initial high RF amplitude
  PulseDur.set_description("Pulse duration of excitation/refocusing pulse");

  EPIRef=false;
  EPIRef.set_description("Interleaved EPI reference scan");

  // register method parameters for user interface, parameter files, etc.
  append_parameter(Blades,"Blades");
  append_parameter(ShortAxis,"ShortAxis");

  append_parameter(TakeMinEchoTime,"TakeMinEchoTime");
  append_parameter(DummyCycles,"DummyCycles");
  append_parameter(RampSampling,"RampSampling");
  append_parameter(FatSaturation,"FatSaturation");
  append_parameter(T1Ernst,"T1Ernst");
  append_parameter(FieldMap,"FieldMap");

  fmapscan.init("fmapscan");
  append_parameter(fmapscan.get_parblock(),"FieldMapPars");


  append_parameter(PulseDur,"PulseDur");

  append_parameter(fMRITrigger,"fMRITrigger");

  append_parameter(UpDownBlade,"UpDownBlade");

  append_parameter(EPIRef,"EPIRef");


  if(systemInfo->get_platform()!=numaris_4) {
    append_parameter(BladeOversampling,"BladeOversampling");
    append_parameter(RampMode,"RampMode");
    append_parameter(RampSteepness,"RampSteepness");
    append_parameter(NumOfGradEchoes,"NumOfGradEchoes",noedit);
    append_parameter(NumOfSamples,"NumOfSamples",noedit);
  }

}



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

  float gamma=systemInfo->get_gamma();

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


  float slicethick=geometryInfo->get_sliceThickness();

  float spatres=slicethick/4.0;

  // calculate Ernst angle accordng to TR
  float flipangle=commonPars->get_FlipAngle();
  if(T1Ernst>0.0) {
    flipangle=180.0/PII * acos( exp ( -secureDivision ( commonPars->get_RepetitionTime(), T1Ernst) ) );
    commonPars->set_FlipAngle( flipangle );
  }

  // excitation pulse
  exc=SeqPulsarSinc("exc", slicethick, true,PulseDur,flipangle,spatres);
  exc.set_rephased(true, 0.8*systemInfo->get_max_grad()); // short rephaser
  exc.set_freqlist( gamma * exc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector() );
  exc.set_pulse_type(excitation);


  // fat saturation module
  fatsat=SeqSat("fatsat",fat);


  //////////////// EPI-Readout: //////////////////////////////

  // square FOV
  int sizeRadial=commonPars->get_MatrixSize(readDirection);
  commonPars->set_MatrixSize(phaseDirection,sizeRadial,noedit);

  commonPars->set_MatrixSize(sliceDirection,1,noedit);


  int readpts_blade=sizeRadial;
  int pelines_blade=sizeRadial;

  int bladewidth=sizeRadial;

  if(Blades>1) {
    float min_bladewidth=sizeRadial*tan(0.5*PII/Blades);  // exact solution
    bladewidth=int( ( (1.0-0.01*BladeOversampling) * min_bladewidth + 0.01*BladeOversampling*sizeRadial ) +0.5 );
  }

  if(ShortAxis) readpts_blade=bladewidth;
  else          pelines_blade=bladewidth;


  float os_read=2.0; // For reduced undersampling artifacts

  float fov=geometryInfo->get_FOV(readDirection); // uniform FOV

  epiacq=SeqAcqEPI("epiacq",commonPars->get_AcqSweepWidth(),
                   readpts_blade, fov,
                   pelines_blade, fov,
                   1, 1, os_read, "", 0, 0, rampType(int(RampMode)),
                   RampSampling, RampSteepness, 0.0);



  // display sampling extents in read/phase direction
  NumOfGradEchoes=epiacq.get_numof_gradechoes();
  NumOfSamples=epiacq.get_npts_read();

  // Template scan fo EPI, begin with copy of actual EPI
  epiacq_template=epiacq;
  epiacq_template.set_label("epiacq_template");

  // 1D phase correction
  epiacq_template.set_template_type(phasecorr_template);


  // Delay instead of actual EPI readout for dummy scans
  epiacq_dummy=SeqDelay("epiacq_dummy",epiacq.get_duration());


  // EPI pre-dephase gradient
  deph=SeqAcqDeph("deph",epiacq);
  deph_template=SeqAcqDeph("deph_template",epiacq_template);




  /////////////////// Rotation of Blades ////////////////////////////////////////////////

  bladerot=SeqRotMatrixVector("bladerot");

  bool full_cycle=UpDownBlade;

  int nangles=Blades;
  if(full_cycle) nangles=2*Blades;

  bladeangels.resize(nangles);
  for(int iangle=0; iangle<nangles; iangle++) {
    RotMatrix rm("rotmatrix"+itos(iangle));
    int iblade=iangle;
    if(full_cycle) iblade=iangle/2;
    float ang=float(iblade)/float(Blades)*PII; // 0...180 deg
    if(full_cycle && iangle%2) ang+=PII; // adjacent blades have opposite direction
    rm.set_inplane_rotation(ang);
    bladeangels[iangle]=ang;
    bladerot.append(rm);
  }


  //////////////// EPI-Reference: //////////////////////////////

  if(EPIRef) {

    // reference gets GRAPPA and partial fourier acquisition
    epiacq_ref=SeqAcqEPI("epiacq_ref",commonPars->get_AcqSweepWidth(),
                   sizeRadial, fov,
                   sizeRadial, fov,
                   1, commonPars->get_ReductionFactor(), os_read, "", 0, 0, rampType(int(RampMode)),
                   RampSampling, RampSteepness, commonPars->get_PartialFourier());


    // Template scan fo EPI, begin with copy of actual EPI
    epiacq_ref_template=epiacq_ref;
    epiacq_ref_template.set_label("epiacq_ref_template");

    // 1D phase correction
    epiacq_ref_template.set_template_type(phasecorr_template);


    // Full multi-shot EPI readout as GRAPPA training data
    epiacq_ref_grappa=epiacq_ref;
    epiacq_ref_grappa.set_label("epiacq_ref_grappa");
    epiacq_ref_grappa.set_template_type(grappa_template);


    // Delay instead of actual EPI readout for dummy scans
    epiacq_ref_dummy=SeqDelay("epiacq_ref_dummy",epiacq_ref.get_duration());


    // EPI pre-dephase gradient
    deph_ref=SeqAcqDeph("deph_ref",epiacq_ref);
    deph_ref_template=SeqAcqDeph("deph_ref_template",epiacq_ref_template);
    deph_ref_grappa=SeqAcqDeph("deph_ref_grappa",epiacq_ref_grappa);
  }


  /////////////////// RF Spoiling ///////////////////////////////////////////////////////

  if(commonPars->get_RFSpoiling()) {

     // recommended by Goerke et al., NMR Biomed. 18, 534-542 (2005)
    int plistsize=16;
    double plistincr=45.0;

    exc.set_phasespoiling(plistsize, plistincr);

    epiacq.set_phasespoiling(plistsize, plistincr);
    epiacq_template.set_phasespoiling(plistsize, plistincr);

    epiacq_ref.set_phasespoiling(plistsize, plistincr);
    epiacq_ref_template.set_phasespoiling(plistsize, plistincr);
    epiacq_ref_grappa.set_phasespoiling(plistsize, plistincr);

    phaseiter=SeqVecIter("phaseiter");

    phaseiter.add_vector(exc.get_phaselist_vector());
    phaseiter.add_vector(epiacq.get_phaselist_vector());
    phaseiter.add_vector(epiacq_template.get_phaselist_vector());

    phaseiter.add_vector(epiacq_ref.get_phaselist_vector());
    phaseiter.add_vector(epiacq_ref_template.get_phaselist_vector());
    phaseiter.add_vector(epiacq_ref_grappa.get_phaselist_vector());

  }

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

  // loop to iterate over slices
  sliceloop=SeqObjLoop("sliceloop");

  // loop to iterate over repetitions
  reploop=SeqObjLoop("reploop");

  // loop to iterate over blades
  bladeloop=SeqObjLoop("bladeloop");

  // loop to iterate over dummy scans
  dummyloop=SeqObjLoop("dummyloop");


  grappaloop_ref=SeqObjLoop("grappaloop_ref");


  //////////////// Timing Delays: //////////////////////////////

  trdelay=SeqDelay("trdelay");

  trdelay_ref=SeqDelay("trdelay_ref");


  //////////////// Crusher Gradient: //////////////////////////////

  double spoiler_strength=0.5*systemInfo->get_max_grad();
  double spoiler_integral=4.0*fabs(deph.get_gradintegral().sum());


  float crusher_integral=2.0*spoiler_integral;
  crusher=SeqGradTrapezParallel("crusher",crusher_integral,crusher_integral,crusher_integral, spoiler_strength);

  crusherdelay=SeqDelay("crusherdelay",0.1); // Small delay to avoid gradient-induced stimulation


  //////////////// trigger: //////////////////////////////

  trigger=SeqTrigger("fmri_trigger",1.0);

  reset=SeqMagnReset("reset");

  //////////////// Field-map template: //////////////////////////////

  if(FieldMap) {
    if(FatSaturation) fmapscan.build_seq(commonPars->get_AcqSweepWidth(),1.0,fatsat); // pass fat saturation on to field-map scan
    else              fmapscan.build_seq(commonPars->get_AcqSweepWidth(),1.0);
  }


  //////////////// Build the sequence: //////////////////////////////

  preppart+=reset;

  // add fat saturation to template and repetitions
  if(FatSaturation) preppart += fatsat;


  dummypart=    preppart + exc + deph          + exc2acq + epiacq_dummy     + crusherdelay + crusher;
  templatepart= preppart + exc + deph_template + exc2acq + epiacq_template  + crusherdelay + crusher;
  imagingpart=  preppart + exc + deph          + exc2acq + epiacq           + crusherdelay + crusher;


  if(EPIRef) {
    templatepart_ref= preppart + exc + deph_ref_template + exc2acq_ref + epiacq_ref_template  + crusherdelay + crusher;
    grappapart_ref=   preppart + exc + deph_ref_grappa   + exc2acq_ref + epiacq_ref_grappa    + crusherdelay + crusher;
    imagingpart_ref=  preppart + exc + deph_ref          + exc2acq_ref + epiacq_ref           + crusherdelay + crusher;
  }


  templatepart.set_gradrotmatrixvector(bladerot);
  imagingpart. set_gradrotmatrixvector(bladerot);


  if(fMRITrigger) slicepart+= trigger; // trigger for PROPELLER interleave

  slicepart_dummy =    sliceloop( dummypart    + trdelay )[exc];
  slicepart_template = sliceloop( templatepart + trdelay )[exc];
  slicepart  +=        sliceloop( imagingpart  + trdelay )[exc];

  if(commonPars->get_RFSpoiling()) {
    slicepart              += phaseiter;
    slicepart_dummy        += phaseiter;
    slicepart_template     += phaseiter;
  }



  if(EPIRef) {

    if(fMRITrigger) slicepart+= trigger; // trigger for EPIRef interleave

    slicepart_template_ref = sliceloop( templatepart_ref + trdelay_ref )[exc];
    slicepart_grappa_ref =   sliceloop( grappapart_ref   + trdelay_ref )[exc];
    slicepart +=             sliceloop( imagingpart_ref  + trdelay_ref )[exc];

    if(commonPars->get_RFSpoiling()) {
      slicepart              += phaseiter;
      slicepart_template_ref += phaseiter;
      slicepart_grappa_ref   += phaseiter;
    }
  }

  if(FieldMap) scan += fmapscan + trdelay;


  if(DummyCycles>0) {
    scan+= dummyloop(
             slicepart_dummy
           )[DummyCycles];
  }


  scan += bladeloop(
            slicepart_template
          )[bladerot];



  if(EPIRef) {

    scan += slicepart_template_ref;

    if(commonPars->get_ReductionFactor()>1) {
      // Fully sampled k-space
      scan+= grappaloop_ref(
               slicepart_grappa_ref
             )[deph_ref_grappa.get_epi_reduction_vector()];
    }
  }


  scan+= reploop( 
           bladeloop(
             slicepart
           )[bladerot]
         )[commonPars->get_NumOfRepetitions()];


  set_sequence( scan );
}




void method_rels() {


   ////////////////// TE Timings: ////////////////////////////////

   double exc_te=exc.get_duration()-exc.get_magnetic_center();

   double min_echo_time_blade= exc_te + deph.get_duration()+epiacq.get_acquisition_center();

   double min_echo_time_ref=0.0;
   if(EPIRef) {
     min_echo_time_ref=exc_te + deph_ref.get_duration()+epiacq_ref.get_acquisition_center();
   }
   double min_echo_time=STD_max(min_echo_time_blade, min_echo_time_ref);

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

   if(TakeMinEchoTime) commonPars->set_EchoTime(min_echo_time);

   exc2acq=    commonPars->get_EchoTime()-min_echo_time_blade;
   exc2acq_ref=commonPars->get_EchoTime()-min_echo_time_ref;




   ////////////////// TR Timings: ////////////////////////////////

   double slicedur    =imagingpart.get_duration()*geometryInfo->get_nSlices();
   double slicedur_ref=imagingpart_ref.get_duration()*geometryInfo->get_nSlices();

   double mintr=STD_max(slicedur, slicedur_ref);

   if(commonPars->get_RepetitionTime()<mintr) commonPars->set_RepetitionTime(mintr);

   trdelay    =secureDivision(commonPars->get_RepetitionTime()-slicedur,     geometryInfo->get_nSlices());
   trdelay_ref=secureDivision(commonPars->get_RepetitionTime()-slicedur_ref, geometryInfo->get_nSlices());

}


void method_pars_set() {

  // extra information for the automatic reconstruction
  epiacq.         set_default_reco_index(userdef,0).set_reco_vector(slice,exc);
  epiacq_template.set_default_reco_index(userdef,0).set_reco_vector(slice,exc);

  epiacq_ref.         set_default_reco_index(userdef,1).set_reco_vector(slice,exc);
  epiacq_ref_template.set_default_reco_index(userdef,1).set_reco_vector(slice,exc);
  epiacq_ref_grappa.  set_default_reco_index(userdef,1).set_reco_vector(slice,exc);

  epiacq.set_reco_vector(cycle,bladerot);
  epiacq_template.set_reco_vector(cycle,bladerot);
  recoInfo->set_DimValues(cycle,bladeangels);

}

};


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


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT