odin 1.8.8-2ubuntu1 / usr / share / odin / sequences / odindti.cpp

#include <odinseq/seqall.h>

class METHOD_CLASS : public SeqMethod {

private:
  JDXfloat  bValue;
  JDXint    NumOfDirections;
  JDXint    b0Rate;
  JDXfloat  GradStrength;

  JDXint    Blades;
  JDXbool   ShortAxis;
  JDXenum   TemplateScan;
  JDXbool   RampSampling;
  JDXenum   RampMode;
  JDXfloat  RampSteepness;
  JDXbool   FatSaturation;
  JDXint    DummyCycles;
  JDXint    NumOfGradEchoes;
  JDXint    NumOfSamples;
  JDXdouble PulseDur;
  JDXbool   FieldMap;


  SeqPulsar  exc;
  SeqPulsar  refoc;
  SeqSat     fatsat;
  SeqAcqEPI  epiacq;
  SeqObjList epipart;
  SeqDelay   epiacq_dummy;
  SeqObjList epipart_dummy;
  SeqAcqEPI  epiacq_template;
  SeqObjList epipart_template;
  SeqAcqEPI  epiacq_grappa;
  SeqObjList epipart_grappa;
  SeqAcqDeph deph;
  SeqAcqDeph deph_template;
  SeqAcqDeph deph_grappa;
  SeqObjLoop sliceloop;
  SeqObjLoop reploop;
  SeqObjLoop bladeloop;
  SeqObjLoop grappaloop;
  SeqObjLoop dummyloop;
  SeqObjLoop dwloop;
  SeqDelay   trdelay;
  SeqObjList scan;
  SeqObjList dummypart;
  SeqObjList templatepart;
  SeqObjList grappapart;
  SeqObjList imagingpart;
  SeqObjList slicepart;
  SeqObjList slicepart_dummy;
  SeqObjList slicepart_template;
  SeqObjList slicepart_grappa;
  SeqObjList preppart;
  SeqDelay   exc2refoc;
  SeqDelay   refoc2acq;
  SeqTrigger trigger;
  SeqGradTrapezParallel spoiler;
  SeqGradTrapezParallel crusher;
  SeqDelay   crusherdelay;
  SeqFieldMap fmapscan;
  SeqVecIter phaseiter;

  SeqDiffWeight dw;
  SeqObjList midpart;

  SeqRotMatrixVector bladerot;

  SeqHalt physiotrigger;
  
  dvector bladeangels;

public:

// This constructor creates an empty EPI sequence
METHOD_CLASS(const STD_string& label) : SeqMethod(label) {
  set_description("Diffusion tensor sequence based on PROPELLER-EPI. A simple Stejskal-Tanner gradient pair around the refocusing pulse is used for diffusion weighting.");
}


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);
  commonPars->set_PhysioTrigger(true);

  bValue=1000.0;
  bValue.set_unit("s/mm^2").set_description("Uniform b-value");

  NumOfDirections=6;
  NumOfDirections.set_description("Number of uniformly distributed diffusion directions");

  b0Rate=5;
  b0Rate.set_description("Rate of b0 scans interleaved wit diffusion weighting");

  GradStrength=80.0;
  GradStrength.set_minmaxval(0.0,100.0).set_unit("%").set_description("Relative strength of diffusion-weighting gradients");

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

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

  TemplateScan.add_item("NoCorrection");
  TemplateScan.add_item("PhaseCorrection");
  TemplateScan.set_actual("PhaseCorrection");
  TemplateScan.set_description("The type of template scan which is acquired beforehand");

  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");

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

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

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


  // register method parameters for user interface, parameter files, etc.
  append_parameter(bValue,"bValue");
  append_parameter(NumOfDirections,"NumOfDirections");
  append_parameter(b0Rate,"b0Rate");
  append_parameter(GradStrength,"GradStrength");

  append_parameter(Blades,"Blades");
  append_parameter(ShortAxis,"ShortAxis");
  append_parameter(DummyCycles,"DummyCycles");
  append_parameter(TemplateScan,"TemplateScan");
  append_parameter(RampSampling,"RampSampling");
  append_parameter(FatSaturation,"FatSaturation");
  append_parameter(FieldMap,"FieldMap");

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

  append_parameter(PulseDur,"PulseDur");

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


}



void method_seq_init() {

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

  float slicethick=geometryInfo->get_sliceThickness();
  float slicegap=geometryInfo->get_sliceDistance()-slicethick;

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

  // Slightly thicker refocusing slice for better SNR
  // Since the same spatial resolution is used for exc and refoc, the gradient strengths will be the same
  float extra_slicethick_refoc=STD_min(0.5*slicethick, 0.3*slicegap);

  // refocusing pulse
  refoc=SeqPulsarSinc("refoc",slicethick+extra_slicethick_refoc,false,PulseDur,180.0);
  refoc.set_freqlist( systemInfo->get_gamma() * refoc.get_strength() / (2.0*PII) * geometryInfo->get_sliceOffsetVector() );
  if(!commonPars->get_RFSpoiling()) refoc.set_phase(90.0);
  refoc.set_pulse_type(refocusing);

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


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

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

  int readpts_blade=sizeRadial;
  int pelines_blade=sizeRadial;

  int bladewidth=sizeRadial;
  if(Blades>1) bladewidth=int(0.5*PII*sizeRadial/Blades+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, commonPars->get_ReductionFactor(), os_read, "", 0, 0, rampType(int(RampMode)),
                   RampSampling, RampSteepness, commonPars->get_PartialFourier());




  // 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
  if(TemplateScan=="PhaseCorrection") {
    epiacq_template.set_template_type(phasecorr_template);
  }


  // Full multi-shot EPI readout as GRAPPA training data
  epiacq_grappa=epiacq;
  epiacq_grappa.set_label("epiacq_grappa");
  epiacq_grappa.set_template_type(grappa_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);
  deph_grappa=SeqAcqDeph("deph_grappa",epiacq_grappa);


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

  bladerot=SeqRotMatrixVector("bladerot");

  bladeangels.resize(Blades);
  for(int iblade=0; iblade<Blades; iblade++) {
    RotMatrix rm("rotmatrix"+itos(iblade));
    float ang=float(iblade)/float(Blades)*PII; // 0...180 deg
    rm.set_inplane_rotation(ang);
    bladeangels[iblade]=ang;
    bladerot.append(rm);
  }



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

  if(commonPars->get_RFSpoiling()) {

     // Use defaults for phase spoiling
    int plistsize=80;
    double plistincr=117.0;

    exc.set_phasespoiling(plistsize, plistincr);
    refoc.set_phasespoiling(plistsize, plistincr, 90.0);

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

    phaseiter=SeqVecIter("phaseiter");
    phaseiter.add_vector(exc.get_phaselist_vector());
    phaseiter.add_vector(refoc.get_phaselist_vector());
    phaseiter.add_vector(epiacq.get_phaselist_vector());
    phaseiter.add_vector(epiacq_template.get_phaselist_vector());
    phaseiter.add_vector(epiacq_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 GRAPPA interleaves to obtain training data
  grappaloop=SeqObjLoop("grappaloop");

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

  // loop to iterate over diffusion weighting
  dwloop=SeqObjLoop("dwloop");

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

  trdelay=SeqDelay("trdelay");


  //////////////// Spoiler Gradient: //////////////////////////////

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

  // on all three axes for uniform b-value
  spoiler=SeqGradTrapezParallel("spoiler",spoiler_integral,spoiler_integral,spoiler_integral, spoiler_strength);


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

  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


  //////////////// 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);
  }

  //////////////// Diffusion weighting: //////////////////////////////

  fvector bvals(1);
  bvals[0]=1000.0*bValue;

  midpart = spoiler + refoc + spoiler;

  dw=SeqDiffWeight("dw", NumOfDirections, bvals, 0.01*GradStrength*systemInfo->get_max_grad(), midpart, b0Rate, true);


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


  // add physiological trigger
  physiotrigger=SeqHalt("physiotrigger");
  if(commonPars->get_PhysioTrigger()) preppart += physiotrigger;


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


  // Rotate deph and epi simultaneously
  epipart_dummy   =SeqObjList("epipart_dummy");
  epipart_template=SeqObjList("epipart_template");
  epipart_grappa  =SeqObjList("epipart_grappa");
  epipart         =SeqObjList("epipart");
  epipart_dummy   += deph          + epiacq_dummy;
  epipart_template+= deph_template + epiacq_template;
  epipart_grappa  += deph_grappa   + epiacq_grappa;
  epipart         += deph          + epiacq;
  epipart_dummy   .set_gradrotmatrixvector(bladerot);
  epipart_template.set_gradrotmatrixvector(bladerot);
  epipart_grappa  .set_gradrotmatrixvector(bladerot);
  epipart         .set_gradrotmatrixvector(bladerot);


  dummypart=    preppart + exc + exc2refoc + dw + refoc2acq + epipart_dummy    + crusherdelay + crusher;

  templatepart= preppart + exc + exc2refoc + dw + refoc2acq + epipart_template + crusherdelay + crusher;

  grappapart=   preppart + exc + exc2refoc + dw + refoc2acq + epipart_grappa   + crusherdelay + crusher;

  imagingpart=  preppart + exc + exc2refoc + dw + refoc2acq + epipart          + crusherdelay + crusher;


  slicepart_dummy =    sliceloop( dummypart + trdelay    )[exc][refoc];

  slicepart_template = sliceloop( templatepart + trdelay )[exc][refoc];

  slicepart_grappa  =  sliceloop( grappapart + trdelay   )[exc][refoc];

  slicepart +=         sliceloop( imagingpart + trdelay  )[exc][refoc];


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

  if(FieldMap) scan += fmapscan + trdelay;


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



  if(TemplateScan=="PhaseCorrection") {
    scan += bladeloop( // Use b0 scan for phase correction of all
              slicepart_template
            )[bladerot];
  }


  if(commonPars->get_ReductionFactor()>1) {
    // Fully sampled k-space
    scan+= grappaloop( // Use b0 scan for interpolation of all
             bladeloop(
               slicepart_grappa
             )[bladerot]
           )[deph_grappa.get_epi_reduction_vector()];
  }


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


  set_sequence( scan );
}




void method_rels() {


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


   double min_echo_time=0.0;



   //////////////// Duration from the middle of the excitation pulse ///////////////////
   //////////////// to the middle of the refocusing pulse: /////////////////////////////

   double TE1=( exc.get_duration() - exc.get_magnetic_center() )
                + dw.get_grad1_duration()
                + spoiler.get_duration()
                + refoc.get_magnetic_center();

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

   //////////////// Duration from the middle of the refocusing pulse ///////////////////
   //////////////// to the middle of the acquisition window: ///////////////////////////

   double TE2=(refoc.get_duration() - refoc.get_magnetic_center())
                + spoiler.get_duration()
                + dw.get_grad2_duration()
                + deph.get_duration()
                + epiacq.get_acquisition_center();

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



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

   commonPars->set_EchoTime(min_echo_time);

   exc2refoc=commonPars->get_EchoTime()/2.0-TE1;
   refoc2acq=commonPars->get_EchoTime()/2.0-TE2;



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

   // reset before recalculating timings
   trdelay=0.0;

   float slicedur=slicepart.get_duration();
   
   if(commonPars->get_PhysioTrigger()) {
     commonPars->set_RepetitionTime(slicedur); // set to minimum TR
   } else {
     if(commonPars->get_RepetitionTime()<slicedur) commonPars->set_RepetitionTime(slicedur);
     trdelay=(commonPars->get_RepetitionTime()-slicedur)/double(geometryInfo->get_nSlices());
   }

}


void method_pars_set() {

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

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

  epiacq.set_reco_vector(dti, dw);
  epiacq_template.set_reco_vector(dti, dw);
  epiacq_grappa.set_reco_vector(dti, dw);
  recoInfo->set_DimValues(dti, dw.get_b_vectors());

}

};


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


// entry point for the sequence module
ODINMETHOD_ENTRY_POINT