Feature module: Dixon (ksdixon.[h,cc])

Data Structures
struct	KSDIXON_DESIGN
struct	KSDIXON_READTRAP
struct	KSDIXON_DUALREADTRAP
struct	KSDIXON_STATE

Macros
#define	KSDIXON_MAX_POINTS   2 /* Maximum number of chemical shift encodes */
#define	KSDIXON_UPDATE_UI(dixon_mode, kdesign, ks_op_dixon_mode, ks_op_dixon_tr_mode, ks_op_dixon_shift)   ksdixon_update_UI(dixon_mode, kdesign, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift))
#define	KSDIXON_INIT_UI(rbw, xres, fov, ks_op_dixon_mode, ks_op_dixon_tr_mode, ks_op_dixon_shift)   ksdixon_init_UI(rbw, xres, fov, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift))
#define	KSDIXON_SETUP(dixon_settings, ks_op_dixon_mode, ks_op_dixon_tr_mode, ks_op_dixon_shift)   ksdixon_setup(dixon_settings, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift));
#define	KSDIXON_MODE_DESCR   "Dixon [1=shift, 2=dBW, 3=spl, 4=invrmp, 5=triptr, 6=dBWs}"
#define	KSDIXON_MODE_NUM_MODES   7
#define	KSDIXON_TR_MODE_DESCR   "TR [0:single 1:multi]"
#define	KSDIXON_SHIFT_DEFAULT   (30000.0/(float)cffield * 1000)
#define	KSDIXON_SHIFT_DESC   "Shift [us]"
#define	KSDIXON_INIT_DESIGN   {KSDIXON_OFF, 1, KS_INITVALUE(KSDIXON_MAX_POINTS,0), KS_NOTSET, 1.0, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KSDIXON_SINGLE_TR}
#define	KSDIXON_INIT_READTRAP   {KS_INITVALUE(4,0.0), KS_INITVALUE(4,0.0)}
#define	KSDIXON_INIT_DUALREADTRAP   {KSDIXON_INIT_READTRAP, KSDIXON_INIT_READTRAP, 0.0, 0.0, 0}
#define	KSDIXON_INIT_STATE   {KS_INITVALUE(3, KS_NOTSET), KS_INITVALUE(KSDIXON_MAX_POINTS,0), 0}

Enumerations
enum	KSDIXON_MODE {   KSDIXON_OFF = 0, KSDIXON_SHIFTED = 1, KSDIXON_DBW = 2, KSDIXON_ASYM_SPLINE = 3,   KSDIXON_ASYM_INVRAMP = 4, KSDIXON_TRIP_TRAP = 5, KSDIXON_DBW_STAGGERED = 6 }
enum	KSDIXON_TR_MODE { KSDIXON_SINGLE_TR = 0, KSDIXON_MULTI_TR = 1 }

Functions
float	ksdixon_calc_dephasing_period ()
float	ksdixon_calc_nsa (float t1, float t2, float var1, float var2, float dephasing_period)
int	ksdixon_max_shift (KSDIXON_MODE mode, int duration, float xres, float fov)
STATUS	ksdixon_dbw_time2echo_hbw (float *time2echo, const float *const t, const float *const G, const float area2echo)
STATUS	ksdixon_dbw_time2echo_lbw (float *time2echo, const float *const t, const float *const G, const float area2echo)
STATUS	ksdixon_time2echo_from_readtrap (float *time2echo, KSDIXON_READTRAP *readtrap, float area2echo)
STATUS	ksdixon_dualreadtrap_eval_shifts (KSDIXON_DUALREADTRAP *dualreadtrap, float area2echo, float shift_at_acq_start, float deadtime)
STATUS	ksdixon_eval_acq_trap (KSDIXON_READTRAP *readtrap, int acq_duration, float sample_area, float max_start_amp, float sampling_threshold_amp, float slew, float amp_limit)
STATUS	ksdixon_eval_dualreadtrap (KSDIXON_DUALREADTRAP *dualreadtrap, int acq_duration1, int acq_duration2, float sample_area, float max_start_end_amp, float sampling_threshold_amp)
STATUS	ksdixon_eval_gre_dBW_acq_waves (KS_WAVE *wave1, KS_WAVE *wave2, int *shift1, int *shift2, int *nover, int acq_start_time, float sampling_threshold_amp, int min_nover, int fov, int res, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_gre_dualBW_readwave (KS_READWAVE *readwave1, KS_READWAVE *readwave2, const KS_KSPACE_DESIGN *kdesign, KSDIXON_DESIGN *dixon, const int acq_start_time, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_asym_spline_acq_wave (KS_WAVE *wave, int shift, int duration, float area, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_trip_trap_acq_points (float *t, float *G, int shift, int duration, float area, float Mc, float g1, float s)
float	ksdixon_dbw_ramptime (const float dur, const float sample_area, const float slewrate)
STATUS	ksdixon_eval_dbw_acq_waves (KS_WAVE *lbw_wave, KS_WAVE *hbw_wave, int *lbw_cse, int *hbw_cse, float *pf, const float min_pf, const int duration, const float area_without_pf, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_flat_acq_wave (KS_WAVE *wave, int duration, float area, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_inverted_ramp_acq_wave (KS_WAVE *wave, int shift, int duration, float area, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_trip_trap_acq_wave (KS_WAVE *wave, int shift, int duration, float area, KS_DESCRIPTION desc)
float	ksdixon_amp_limit (float board_amp_limit, float fov)
STATUS	ksdixon_eval_inverted_ramps (KS_WAVE *acq_waves, int *shifts, int points, const KS_KSPACE_DESIGN *kdesign, KS_DESCRIPTION desc)
STATUS	ksdixon_eval_inverted_ramp_acq_points (float *t, float *G, int shift, int duration, float area, float s, float amp_limit)
void	ksdixon_init_settings (KSDIXON_DESIGN *const dixon)
STATUS	ksdixon_setup (KSDIXON_DESIGN *pdixon, const int opnum_dixon_mode, _cvfloat *_opuser_dixon_mode, const int opnum_dixon_tr_mode, _cvfloat *_opuser_dixon_tr_mode, const int opnum_dixon_shift, _cvfloat *_opuser_dixon_shift)
STATUS	ksdixon_eval_validate_settings (const KSDIXON_DESIGN *dixon)
void	ksdixon_update_UI (const KSDIXON_MODE dixon_mode, const KS_KSPACE_DESIGN *const kdesign, const int opnum_dixon_mode, _cvfloat *_opuser_dixon_mode, const int opnum_dixon_tr_mode, _cvfloat *_opuser_dixon_tr_mode, const int opnum_dixon_shift, _cvfloat *_opuser_dixon_shift)
void	ksdixon_init_UI (float rbw, int xres, float fov, const int opnum_dixon_mode, _cvfloat *_opuser_dixon_mode, const int opnum_dixon_tr_mode, _cvfloat *_opuser_dixon_tr_mode, const int opnum_dixon_shift, _cvfloat *_opuser_dixon_shift)
STATUS	ksdixon_eval_setup_gre_readouts (KS_ECHOTRAIN *echotrain, KSDIXON_DESIGN *dixon, const KS_KSPACE_DESIGN *kdesign, const int acq_start_time, const int TR_index, KS_DESCRIPTION desc)
STATUS	ksdixon_set_design (KSDIXON_DESIGN *const dixon_design, const KS_KSPACE_DESIGN *const kdesign)
STATUS	ksdixon_setup_fse_echotrain (KS_ECHOTRAIN_DESIGN *const echotrain_design, KSDIXON_STATE *const dixon_state, const KSDIXON_DESIGN *const dixon_design, const int ETL, float xcrusher_area, const KS_DESCRIPTION description)
STATUS	ksdixon_eval_trip_traps (KS_WAVE *acq_waves, int *shifts, int points, const KS_KSPACE_DESIGN *kdesign, KS_DESCRIPTION desc)

Variables
int	piuset

Detailed Description

Macro Definition Documentation

KSDIXON_MAX_POINTS

#define KSDIXON_MAX_POINTS   2 /* Maximum number of chemical shift encodes */

KSDIXON_UPDATE_UI

#define KSDIXON_UPDATE_UI	(		dixon_mode,
			kdesign,
			ks_op_dixon_mode,
			ks_op_dixon_tr_mode,
			ks_op_dixon_shift
	)		ksdixon_update_UI(dixon_mode, kdesign, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift))

Macro that calls ksdixon_update_UI. Currently not used. The idea is to update dixon opusers according to kdesign and dixon_mode.

KSDIXON_INIT_UI

#define KSDIXON_INIT_UI	(		rbw,
			xres,
			fov,
			ks_op_dixon_mode,
			ks_op_dixon_tr_mode,
			ks_op_dixon_shift
	)		ksdixon_init_UI(rbw, xres, fov, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift))

Macro that calls ksdixon_init_UI. The purpose of these macros is to directly pass the opuser number rather than the full cv (the macro does the expansion for you). For instance:

In implementation e file

// Select opuser numbers for dixon options
@global
  #define OPUSER_NUM_DIXON_MODE 17
  #define OPUSER_NUM_DIXON_TR_MODE 18
  #define OPUSER_NUM_DIXON_SHIFT 19

ksfse2_init_UI():
  // Setup the opuser description, range, etc.
  KSDIXON_INIT_UI(oprbw, opxres, opfov, OPUSER_NUM_DIXON_MODE, OPUSER_NUM_DIXON_TR_MODE, OPUSER_NUM_DIXON_SHIFT);

Parameters

[in]	rbw	Readout bandwidth (oprbw)
[in]	xres	Number of readout samples prescribed (opxres)
[in]	fov	Readout FOV (opfov) [mm]
[in]	ks_op_dixon_mode	opuser number for selecting dixon mode (e.g. 17 for opuser17)
[in]	ks_op_dixon_tr_mode	opuser number for selecting dixon TR mode (multi/single TR)
[in]	ks_op_dixon_shift	opuser number for selecting dixon cse (shift in us)

KSDIXON_SETUP

#define KSDIXON_SETUP	(		dixon_settings,
			ks_op_dixon_mode,
			ks_op_dixon_tr_mode,
			ks_op_dixon_shift
	)		ksdixon_setup(dixon_settings, KS_EXPAND_OPUSER(ks_op_dixon_mode), KS_EXPAND_OPUSER(ks_op_dixon_tr_mode), KS_EXPAND_OPUSER(ks_op_dixon_shift));

Macro that calls ksdixon_setup. The purpose of these macros is to directly pass the opuser number rather than the full cv (the macro does the expansion for you). For instance:

In implementation e file

// Select opuser numbers for dixon options
@global
  #define OPUSER_NUM_DIXON_MODE 17
  #define OPUSER_NUM_DIXON_TR_MODE 18
  #define OPUSER_NUM_DIXON_SHIFT 19

ksfse2_setup_stretched():
  // Fill ksfse2_dixon_settings with dixon related values.
  status = KSDIXON_SETUP(&ksfse2_dixon_settings, OPUSER_NUM_DIXON_MODE, OPUSER_NUM_DIXON_TR_MODE, OPUSER_NUM_DIXON_SHIFT);

Parameters

[out]	dixon_settings	- An instance of KSDIXON_DESIGN that will be filled according to the values in the opusers
[in]	ks_op_dixon_mode	- opuser number for selecting dixon mode (e.g. 17 for opuser17)
[in]	ks_op_dixon_tr_mode	- opuser number for selecting dixon TR mode (multi/single TR)
[in]	ks_op_dixon_shift	- opuser number for selecting dixon cse (shift in us)

Return values

STATUS

SUCCESS or FAILURE

KSDIXON_MODE_DESCR

#define KSDIXON_MODE_DESCR   "Dixon [1=shift, 2=dBW, 3=spl, 4=invrmp, 5=triptr, 6=dBWs}"

KSDIXON_MODE_NUM_MODES

#define KSDIXON_MODE_NUM_MODES   7

KSDIXON_TR_MODE_DESCR

#define KSDIXON_TR_MODE_DESCR   "TR [0:single 1:multi]"

KSDIXON_SHIFT_DEFAULT

#define KSDIXON_SHIFT_DEFAULT   (30000.0/(float)cffield * 1000)

KSDIXON_SHIFT_DESC

#define KSDIXON_SHIFT_DESC   "Shift [us]"

KSDIXON_INIT_DESIGN

#define KSDIXON_INIT_DESIGN   {KSDIXON_OFF, 1, KS_INITVALUE(KSDIXON_MAX_POINTS,0), KS_NOTSET, 1.0, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KS_NOTSET, KSDIXON_SINGLE_TR}

KSDIXON_INIT_READTRAP

#define KSDIXON_INIT_READTRAP   {KS_INITVALUE(4,0.0), KS_INITVALUE(4,0.0)}

KSDIXON_INIT_DUALREADTRAP

#define KSDIXON_INIT_DUALREADTRAP   {KSDIXON_INIT_READTRAP, KSDIXON_INIT_READTRAP, 0.0, 0.0, 0}

KSDIXON_INIT_STATE

#define KSDIXON_INIT_STATE   {KS_INITVALUE(3, KS_NOTSET), KS_INITVALUE(KSDIXON_MAX_POINTS,0), 0}

Enumeration Type Documentation

KSDIXON_MODE

enum KSDIXON_MODE

Enums for selecting a Dixon mode. This typically changes the readout window but is not restricted to only that.

• KSDIXON_OFF - No Dixon
• KSDIXON_SHIFTED - Conventional CSE using shifted gradients [Hardy PA, Hinks RS, Tkach JA. Separation of fat and water in fast spin‐echo MR imaging with the three‐point dixon technique. J. Magn. Reson. Imaging 1995;5:181–185.]
• KSDIXON_DBW - Dual bandwidths [Rydén H, Berglund J, Norbeck O, et al. RARE two-point Dixon with dual bandwidths. Magn. Reson. Med. 2020;84:2456–2468.]
• KSDIXON_ASYM_SPLINE - Asymmetric spline [Rydén H, Norbeck O, Avventi E, et al. Chemical shift encoding using asymmetric readout waveforms. Magn. Reson. Med. 2020;42:963.]
• KSDIXON_ASYM_INVRAMP - Maximum CSE
• KSDIXON_TRIP_TRAP - Asymmetric piece-wise linear waveform with three bandwidths [FSE Dixon with bandwidth-matched asymmetric readouts tailored for real-valued fat/water estimates, ISMRM 2023]

Enumerator
KSDIXON_OFF
KSDIXON_SHIFTED
KSDIXON_DBW
KSDIXON_ASYM_SPLINE
KSDIXON_ASYM_INVRAMP
KSDIXON_TRIP_TRAP
KSDIXON_DBW_STAGGERED

             {KSDIXON_OFF = 0,
              KSDIXON_SHIFTED = 1,
              KSDIXON_DBW = 2,
              KSDIXON_ASYM_SPLINE = 3,
              KSDIXON_ASYM_INVRAMP = 4,
              KSDIXON_TRIP_TRAP = 5,
              KSDIXON_DBW_STAGGERED = 6
} KSDIXON_MODE;

KSDIXON_TR_MODE

enum KSDIXON_TR_MODE

Enumerator
KSDIXON_SINGLE_TR
KSDIXON_MULTI_TR

             {KSDIXON_SINGLE_TR = 0, /* Echoes in echo train toggle between chemical shift encodes */
              KSDIXON_MULTI_TR = 1 /* One CSE per shot */
} KSDIXON_TR_MODE;

Function Documentation

ksdixon_calc_dephasing_period()

float ksdixon_calc_dephasing_period

(

)

Period of dephasing between fat and water [usec]

Uses cv cffield to determine field strength.

Returns

[us] float

                                      {
  float delta = 3.4e-6; /* chemical shift is 3.4 ppm */
  float gyro = 42.5759; /* MHz/T */
  float B0 = cffield * 1e-4; /* T */

  return 1.0f / (gyro * B0 * delta); /* usec */

} /* ksdixon_calc_dephasing_period() */

ksdixon_calc_nsa()

float ksdixon_calc_nsa	(	float	t1,
		float	t2,
		float	var1,
		float	var2,
		float	dephasing_period
	)

Calculate mean single-species weighted NSA for 2-point Dixon

Following Eq. 7 in Ryden et al. MRM 84(5):2456-68, 2020 (using noise variance var=1/w^2)

Parameters

t1
t2
var1
var2
dephasing_period

Returns

float

                                                                                           {
    float wsq1 = 1.0 / var1;
    float wsq2 = 1.0 / var2;
    const float c1 = cos(2 * PI / dephasing_period * t1);
    const float c2 = cos(2 * PI / dephasing_period * t2);

    return (wsq1 * wsq2 * pow(c1 - c2, 2) * (wsq1*c1*c1 + wsq2*c2*c2 + wsq1 + wsq2) / (2 * (wsq1 + wsq2) * (wsq1*c1*c1 + wsq2*c2*c2)) );

} /* ksdixon_calc_nsa() */

ksdixon_max_shift()

int ksdixon_max_shift	(	KSDIXON_MODE	mode,
		int	duration,
		float	xres,
		float	fov
	)

Find the max shift for for a certain Dixon mode using binary sarch. Used for get ranges in the UI

Parameters

[in]	mode	The dixon mode
[in]	duration	[us] Acq duration
[in]	xres	Number of samples along x
[in]	fov	[mm] FOV

Returns

int [us] max CSE

                                                                              {
  STATUS status;
  float sample_area = xres * ks_calc_fov2gradareapixel(fov);
  float tsp = duration / xres;
  float trap_amp = ((float) (1.0 / (GAM * tsp * 1e-6) * (10.0 / fov)));

  KS_WAVE acq_wave;
  int right_shift = duration / 2; /* Maximum shift, shouldn't be possible */
  int left_shift = 0;
  int shift;
  int iters = 20;
  int diff = duration;
  int i = 0;
  /* Binary search */
  for (; (i < iters) && (diff >= 8) /* us */; i++) {
    shift = (right_shift + left_shift) / 2;
    if (mode == KSDIXON_ASYM_SPLINE) {
      KS_DESCRIPTION tmpdesc = "ui_maxshift";
      status = ksdixon_eval_asym_spline_acq_wave(&acq_wave, shift, duration, sample_area, tmpdesc);
    } else if (mode == KSDIXON_TRIP_TRAP) {
      float t[6] = {KS_NOTSET};
      float G[6] = {KS_NOTSET};
      status = ksdixon_eval_trip_trap_acq_points(t, G, shift, duration, sample_area, .25*sample_area, trap_amp, ks_syslimits_slewrate1(loggrd));
    } else if (mode == KSDIXON_ASYM_INVRAMP) {
      float t[4] = {};
      float G[4] = {};
      status = ksdixon_eval_inverted_ramp_acq_points(t, G, shift, duration, sample_area, .95 * ks_syslimits_slewrate1(loggrd), ksdixon_amp_limit(ks_syslimits_ampmax1(loggrd), fov));
    } else if (mode == KSDIXON_DBW || mode == KSDIXON_DBW_STAGGERED) {
      return 0;
    } else {
      return ks_error("%s: Not implemented for the requested dixon method. It is easy and you should do it now.", __FUNCTION__);
    }
    if (status == SUCCESS) {
      diff = shift - left_shift;
      left_shift = shift;
    } else {
      right_shift = shift;
    }
  }
  return left_shift;

} /* ksdixon_max_shift() */

ksdixon_dbw_time2echo_hbw()

STATUS ksdixon_dbw_time2echo_hbw	(	float *	time2echo,
		const float *const	t,
		const float *const	G,
		const float	area2echo
	)

                                                                                                                        {
  const float ramp_up_area = (t[2] - t[0]) * (G[2] + G[0]) / 2.0;
  const float plateau_area = (t[3] - t[2]) * G[2];
  if (area2echo < ramp_up_area) { /* echo is on ramp up */
    const float slew = (G[2] - G[0]) / (t[2] - t[0]);
    *time2echo = sqrt(2 * area2echo / slew);
  } else if (area2echo < (ramp_up_area + plateau_area) ) { /* echo is on plateau */ 
    *time2echo = t[2] + (area2echo - ramp_up_area) / G[2];
  } else {
    return ks_error("%s: hbw cannot reach area2echo (%f)", __FUNCTION__, area2echo);
  }

  return SUCCESS;

} /* ksdixon_dbw_time2echo_hbw() */

ksdixon_dbw_time2echo_lbw()

STATUS ksdixon_dbw_time2echo_lbw	(	float *	time2echo,
		const float *const	t,
		const float *const	G,
		const float	area2echo
	)

                                                                                                                        {
  const float ramp_up_area = (t[1] - t[0]) * (G[1] + G[0]) / 2.0;
  const float plateau_area = (t[2] - t[1]) * G[1];
  if (area2echo < (ramp_up_area + plateau_area)) { /* echo is on plateau */
    *time2echo = t[1] + (area2echo - ramp_up_area) / G[1];
  } else { /* echo is on ramp down */ 
    const float slew = (G[2] - G[3]) / (t[3] - t[2]);
    const float square = pow((G[2] / slew), 2) - 2.0 * (area2echo - ramp_up_area - plateau_area) / slew;
    if (square < 0.0) {
      return ks_error("%s: lbw cannot reach area2echo (%f)", __FUNCTION__, area2echo);
    }
    *time2echo = t[2] + G[2] / slew - sqrt(square);
  }

  return SUCCESS;

} /* ksdixon_dbw_time2echo_lbw() */

ksdixon_time2echo_from_readtrap()

STATUS ksdixon_time2echo_from_readtrap	(	float *	time2echo,
		KSDIXON_READTRAP *	readtrap,
		float	area2echo
	)

Calculate the time until k-space center is reached given a KSDIXON_READTRAP and the area to center

Parameters

[out]	time2echo	[us]
[in]	readtrap	A readtrap support point struct
[in]	area2echo	Area required to reach k-space center

Return values

STATUS

SUCCESS or FAILURE

                                                                                                      {
  float* t = readtrap->t;
  float* G = readtrap->G;

  return ksdixon_dbw_time2echo_lbw(time2echo, t, G, area2echo);

} /* ksdixon_time2echo_from_readtrap() */

ksdixon_dualreadtrap_eval_shifts()

STATUS ksdixon_dualreadtrap_eval_shifts	(	KSDIXON_DUALREADTRAP *	dualreadtrap,
		float	area2echo,
		float	shift_at_acq_start,
		float	deadtime
	)

Calculate the shifts in a Dixon dbw readout.

Parameters

[out]	dualreadtrap	Will set dualreadtrap->shift1 and dualreadtrap->shift2
[in]	area2echo	Area to k-space center
[in]	shift_at_acq_start	[us] CSE when the acq starts
[in]	deadtime	[us] Deadtime corresponds to inner ramps where amp is below sampling threshold

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                       {
  STATUS status;
  float time2echo;

  status = ksdixon_time2echo_from_readtrap(&time2echo, &dualreadtrap->trap1, area2echo);
  KS_RAISE(status);
  dualreadtrap->shift1 = shift_at_acq_start + time2echo;

  status = ksdixon_time2echo_from_readtrap(&time2echo, &dualreadtrap->trap2, area2echo);
  KS_RAISE(status);
  /* note that trap2 is "mirrored" */
  dualreadtrap->shift2 = shift_at_acq_start + dualreadtrap->trap1.t[3] + deadtime + dualreadtrap->trap2.t[3] - time2echo;

  return SUCCESS;

} /* ksdixon_dualreadtrap_eval_shifts() */

ksdixon_eval_acq_trap()

STATUS ksdixon_eval_acq_trap	(	KSDIXON_READTRAP *	readtrap,
		int	acq_duration,
		float	sample_area,
		float	max_start_amp,
		float	sampling_threshold_amp,
		float	slew,
		float	amp_limit
	)

Evaluates readout support points for a trapezoid

Parameters

[out]	readtrap	Struct with support points
[in]	acq_duration	[us] Duration of acquisition window
[in]	sample_area	Sample area
[in]	max_start_amp	Maximum amp at the start of acquisition
[in]	sampling_threshold_amp	Amplitude threshold where sampling is toggled
[in]	slew	System slew rate
[in]	amp_limit	System gradient amplitude limit

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                                              {
  float amp; /* plateau amplitude */
  float rem_area = sample_area; /* remaining sample area */
  rem_area -= acq_duration * sampling_threshold_amp; /* base area up to sampling threshold amp */
  if (rem_area < 0.0) {
    /* ks_dbg("%s: Trapezoid did not reach the sampling threshold", __FUNCTION__); */
    return FAILURE;
  }
  /* "middle area" above base area up to max start amp includes a ramp down */
  float middle_area = (max_start_amp - sampling_threshold_amp) * (acq_duration - (max_start_amp - sampling_threshold_amp) / (2.0 * slew));
  if (rem_area < middle_area) { /* Need to lower the amplitude */
    float square = pow(acq_duration * slew, 2) - 2.0 * rem_area * slew;
    if (square < 0.0) {
      /* ks_dbg("%s: Trapezoid did not reach max_start_amp", __FUNCTION__); */
      return FAILURE;
    }
    amp = acq_duration * slew - sqrt(square) + sampling_threshold_amp;
  } else {
    rem_area -= middle_area; /* remaining "top area" above base and middle areas includes a ramp up and a ramp down */
    float term = (acq_duration * slew - max_start_amp + sampling_threshold_amp) / 2.0;
    float square = pow(term, 2) - rem_area * slew;
    if (square < 0.0) {
      /* ks_dbg("%s: Trapezoid did not reach plateau", __FUNCTION__); */
      return FAILURE;
    }
    amp = term - sqrt(square) + max_start_amp;
  }

  if (amp > amp_limit) {
    /* ks_dbg("%s: Trapezoid amplitude (%f) exceeds amp_limit (%f)", __FUNCTION__, amp, amp_limit); */
    return FAILURE;
  }

  readtrap->t[0] = 0.0;
  readtrap->G[0] = FMin(2, max_start_amp, amp);

  readtrap->t[1] = (amp - readtrap->G[0]) / slew;
  readtrap->G[1] = amp;

  readtrap->t[2] = acq_duration - (amp - sampling_threshold_amp) / slew;
  readtrap->G[2] = amp;

  readtrap->t[3] = acq_duration;
  readtrap->G[3] = sampling_threshold_amp;

  return SUCCESS;

} /* ksdixon_eval_acq_trap() */

ksdixon_eval_dualreadtrap()

STATUS ksdixon_eval_dualreadtrap	(	KSDIXON_DUALREADTRAP *	dualreadtrap,
		int	acq_duration1,
		int	acq_duration2,
		float	sample_area,
		float	max_start_end_amp,
		float	sampling_threshold_amp
	)

Evaluates a dbw readout (Only evaluation, this should be wrapped around some optimization)

Parameters

[out]	dualreadtrap
[in]	acq_duration1	[us] Duration of first acquisition
[in]	acq_duration2	[us] Duration of second acquisition
[in]	sample_area	Sample area for one readout. This value must take partial Fourier into account
[in]	max_start_end_amp	Maximum amp at the start and end of acquisition
[in]	sampling_threshold_amp	Amplitude threshold where sampling is toggled on/off

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                                                     {
  STATUS status;

  float slew_limit = 0.98 * ks_syslimits_slewrate1(loggrd);
  float amp_limit = ks_syslimits_ampmax1(loggrd);

  if (max_start_end_amp > amp_limit) {
    return ks_error("%s: max_start_end_amp (%f) exceeds amp_limit (%f)", __FUNCTION__, max_start_end_amp, amp_limit);
  } 
  if (sampling_threshold_amp > max_start_end_amp) {
    return ks_error("%s: max_start_end_amp (%f) is below sampling_threshold_amp (%f)", __FUNCTION__, max_start_end_amp, sampling_threshold_amp);
  }

  status = ksdixon_eval_acq_trap(&dualreadtrap->trap1, acq_duration1, sample_area, max_start_end_amp, sampling_threshold_amp, slew_limit, amp_limit);
  if (status != SUCCESS) return status;

  status = ksdixon_eval_acq_trap(&dualreadtrap->trap2, acq_duration2, sample_area, sampling_threshold_amp, sampling_threshold_amp, slew_limit, amp_limit);
  if (status != SUCCESS) return status;

  return SUCCESS;

} /* ksdixon_eval_dualreadtrap() */

ksdixon_eval_gre_dBW_acq_waves()

STATUS ksdixon_eval_gre_dBW_acq_waves	(	KS_WAVE *	wave1,
		KS_WAVE *	wave2,
		int *	shift1,
		int *	shift2,
		int *	nover,
		int	acq_start_time,
		float	sampling_threshold_amp,
		int	min_nover,
		int	fov,
		int	res,
		KS_DESCRIPTION	desc
	)

Generates two acq waves for GRE dBW using a brute-force search

Parameters

[out]	wave1	Output wave 1
[out]	wave2	Output wave 2
[out]	shift1	CSE for wave 1
[out]	shift2	CSE for wave 2
[out]	nover	Partial Fourier
[in]	acq_start_time	Time when first acquisition is turned on
[in]	sampling_threshold_amp	Amplitude threshold where sampling is toggled on/off
[in]	min_nover	Minimum number of samples after k-space center ( partial Fourier )
[in]	fov	[mm]
[in]	res	Number of samples in readout direction (as if no partial Fourier)
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                                                                                    {
  STATUS status = FAILURE;

  KSDIXON_DUALREADTRAP dualreadtrap = KSDIXON_INIT_DUALREADTRAP;
  KSDIXON_DUALREADTRAP best_dualreadtrap = KSDIXON_INIT_DUALREADTRAP;

  /* Brute-force optimize with free parameters nover (partial Fourier in read dir) and acq_duration of first trap */
  float max_nsa = 0.0f; /* target for optimization */

  float dephasing_period = ksdixon_calc_dephasing_period(); /* usec */

  /* minimum acq_duration of first trap corresponds to in-phase echo for min_nover */
  int min_dur = RUP_FACTOR((int) ks_round(dephasing_period) / 4, 16 * GRAD_UPDATE_TIME);
  int max_dur = RUP_GRD((int) ks_round(dephasing_period));
  int acq_duration1, acq_duration2;
  for (acq_duration1 = min_dur; acq_duration1 < max_dur; acq_duration1 += 16 * GRAD_UPDATE_TIME) {

    for (acq_duration2 = min_dur; acq_duration2 < max_dur; acq_duration2 += 16 * GRAD_UPDATE_TIME) {

      for (dualreadtrap.nover = min_nover; dualreadtrap.nover <= (res/2); dualreadtrap.nover += 4) {

        float sample_area = (res/2 + dualreadtrap.nover) * ks_calc_fov2gradareapixel(fov);

        if (SUCCESS != ksdixon_eval_dualreadtrap(&dualreadtrap, acq_duration1, acq_duration2, sample_area, sampling_threshold_amp, sampling_threshold_amp)) {
          continue;
        }

        float area2echo = (res/2) * ks_calc_fov2gradareapixel(fov);
        /* deadtime corresponds to inner ramps where amp is below sampling threshold */
        int deadtime = 2 * RUP_GRD((int) (sampling_threshold_amp / ks_syslimits_slewrate1(loggrd)));
        if (SUCCESS != ksdixon_dualreadtrap_eval_shifts(&dualreadtrap, area2echo, (float) acq_start_time, (float) deadtime)) {
          return ks_error("%s: failed to eval shifts of dualreadtrap", __FUNCTION__);
        }

        if (fabs(dualreadtrap.shift1 - dualreadtrap.shift2) > dephasing_period) {
          break; /* if this failed, higher nover will also fail */
        }

        float nsa = ksdixon_calc_nsa(dualreadtrap.shift1, dualreadtrap.shift2, 1.0 / (float) acq_duration1, 1.0 / (float) acq_duration2, dephasing_period);

        if (nsa > max_nsa) {
          max_nsa = nsa;
          best_dualreadtrap = dualreadtrap;
          status = SUCCESS;
        }
      }
    }
  }

  if (status != SUCCESS) {
    return ks_error("%s: unable to find dualreadtrap", __FUNCTION__);
  }

  *shift1 = ks_round(best_dualreadtrap.shift1);
  *shift2 = ks_round(best_dualreadtrap.shift2);
  *nover = best_dualreadtrap.nover;

  KS_DESCRIPTION tmpdesc;
  sprintf(tmpdesc, "%s_readtrap_1", desc);
  status = ks_eval_coords2wave(wave1, best_dualreadtrap.trap1.t, best_dualreadtrap.trap1.G, 4, GRAD_UPDATE_TIME, tmpdesc);
  KS_RAISE(status);

  sprintf(tmpdesc, "%s_readtrap_2", desc);
  status = ks_eval_coords2wave(wave2, best_dualreadtrap.trap2.t, best_dualreadtrap.trap2.G, 4, GRAD_UPDATE_TIME, tmpdesc);
  KS_RAISE(status);

  status = ks_eval_mirrorwave(wave2);
  KS_RAISE(status);

  return SUCCESS;

} /* ksdixon_eval_gre_dBW_acq_waves() */

ksdixon_eval_gre_dualBW_readwave()

STATUS ksdixon_eval_gre_dualBW_readwave	(	KS_READWAVE *	readwave1,
		KS_READWAVE *	readwave2,
		const KS_KSPACE_DESIGN *	kdesign,
		KSDIXON_DESIGN *	dixon,
		const int	acq_start_time,
		KS_DESCRIPTION	desc
	)

Convenience wrapper that calls ksdixon_eval_gre_dBW_acq_waves and generates readwaves

Parameters

[out]	readwave1	Output readwave 1
[out]	readwave2	Output readwave 2
[in]	kdesign	A KS_KSPACE_DESIGN, required to calculate readout areas and nover
[out]	dixon	A KSDIXON_DESIGN object that gets filled with shifts TODO
	acq_start_time	Time when first acquisition is turned on
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                                                               {
  STATUS status;

  int nover;
  KS_WAVE acq_wave1, acq_wave2;

  /* use amplitude of the default trap as sampling threshold (i.e. correspondning to the kdesign.rbw) */
  /* this preserves the acq_delay, as not to corrupt calculations of shifts */
  KS_READTRAP default_trap;
  status = ks_eval_design_readtrap(&default_trap, kdesign, 0.0, "default");
  KS_RAISE(status);
  float sampling_threshold_amp = default_trap.grad.amp;

  int min_nover = (kdesign->nover[XGRAD] > 0) ? kdesign->nover[XGRAD]
                                              : 0;

  status = ksdixon_eval_gre_dBW_acq_waves(&acq_wave1, &acq_wave2, &dixon->shifts[0], &dixon->shifts[1], &nover, acq_start_time, sampling_threshold_amp, min_nover, kdesign->fov[XGRAD], kdesign->res[XGRAD], desc);
  KS_RAISE(status);

  readwave1->fov = kdesign->fov[XGRAD];
  readwave1->res = kdesign->res[XGRAD];
  readwave1->nover = -nover; /* first trap is a "late echo" */

  status = ks_eval_readwave(readwave1, &acq_wave1, 0.0, 0.0, 0, KS_CRUSHER_STRATEGY_FASTEST, KS_CRUSHER_STRATEGY_FASTEST);
  KS_RAISE(status);

  readwave2->fov = kdesign->fov[XGRAD];
  readwave2->res = kdesign->res[XGRAD];
  readwave2->nover = nover; /* second trap is an "early echo" */

  status = ks_eval_readwave(readwave2, &acq_wave2, 0.0, 0.0, 0, KS_CRUSHER_STRATEGY_FASTEST, KS_CRUSHER_STRATEGY_FASTEST);
  KS_RAISE(status);  

  return SUCCESS;

} /* ksdixon_eval_gre_dualBW_readwave() */

ksdixon_eval_asym_spline_acq_wave()

STATUS ksdixon_eval_asym_spline_acq_wave	(	KS_WAVE *	wave,
		int	shift,
		int	duration,
		float	area,
		KS_DESCRIPTION	desc
	)

Evaluate an asymmetric spline-based acq wave according to [Rydén H, et al. Chemical shift encoding using asymmetric readout waveforms. Magn. Reson. Med. 2020;42:963]

Assumes there is not other waveform played (calls ks_syslimits_ampmax1(loggrd))

Parameters

[out]	wave	Output wave
[in]	shift	[us] Desired CSE
[in]	duration	Duration of the acq wave (i.e. sampling duration)
[in]	area	Sample area
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                  {
  STATUS status;

  if (shift > 0) {
    status = ksdixon_eval_asym_spline_acq_wave(wave, -shift, duration, area, desc);
    KS_RAISE(status);
    return ks_eval_mirrorwave(wave);
  } else if (shift == 0) { /* Make a flat acq wave if in-phase echo is desired */
    float t[2] = {0, (float)duration};
    float G[2] = {area/duration, area/duration};
    status = ks_eval_coords2wave(wave, t, G, 2, GRAD_UPDATE_TIME, desc);
    return status;
  }

  if (duration % GRAD_UPDATE_TIME) {
    return ks_error("%s: duration (%d) must fit on the gradient raster (%d)", __FUNCTION__, duration, GRAD_UPDATE_TIME);
  }

  float base_area = area / 2; /* put half the sample area in a base rectangle */

  const int ttc = RDN_GRD(duration / 2 + shift); /* time to center (i.e. gradient echo) */

  const double g0 = base_area / duration;
  const double T1 = ttc;
  const double T2 = duration - T1;
  const double M_b1 = T1/(T1+T2) * base_area;
  const double M_b2 = base_area - M_b1;
  const double M1 = area / 2 - M_b1;
  const double M2 = area / 2 - M_b2;

  const double coeffs_a[4] = { /* a3 */  4*(  M1*T2*T2*T2*T2   - 2*M2*T1*T1*T1*T1   + 6*M1*T1*T1*T2*T2 - 4*M2*T1*T1*T2*T2 +  5*M1*T1*T2*T2*T2 - 8*M2*T1*T1*T1*T2) / (  T1*T1*T1*T1*T2    * (T1 + T2)     * (3*T2*T2 + 2*T1*T2)),
                               /* a2 */ 12*(  M2*T1*T1*T1      + 4*M2*T1*T1*T2      - 4*M1*T1*T2*T2    - M1*T2*T2*T2)                                             / (2*T1*T1*T1*T1*T2*T2 + 3*T1*T1*T1*T2*T2*T2),
                               /* a1 */  4*(3*M1*T2*T2*T2*T2   -   M2*T1*T1*T1*T1   + 6*M1*T1*T1*T2*T2 - 6*M2*T1*T1*T2*T2 + 10*M1*T1*T2*T2*T2 - 6*M2*T1*T1*T1*T2) / (        T1*T1*T2*T2 * (2*T1 + 3*T2) * (T1+T2) ),
                               /* a0 */  g0
                              };
  const double coeffs_b[4] = { -(3*coeffs_a[0]*T1       +   coeffs_a[1]) / (3*T2),                  /* b3 in paper */
                                (3*coeffs_a[0]*T1       +   coeffs_a[1]),                         /* b2 */
                                 3*coeffs_a[0]*T1*T1    + 2*coeffs_a[1]*T1    + coeffs_a[2],      /* b1 */
                                   coeffs_a[0]*T1*T1*T1 +   coeffs_a[1]*T1*T1 + coeffs_a[2]*T1 + coeffs_a[3] }; /* b0 */

  const int num_points_a = ttc / GRAD_UPDATE_TIME;
  const int num_points_b = (duration - ttc) / GRAD_UPDATE_TIME;
  float ta[num_points_a];
  float tb[num_points_b];
  int idx;
  for (idx = 0; idx < num_points_a; idx++) {
    ta[idx] = (idx + 0.5) * GRAD_UPDATE_TIME;
  }
  for (idx = 0; idx < num_points_b; idx++) {
    tb[idx] = (idx + 0.5) * GRAD_UPDATE_TIME;
  }

  KS_WAVEFORM waveform = KS_INIT_WAVEFORM;

  ks_polyval_f(coeffs_a, 3, ta, num_points_a, &waveform[0]);
  ks_polyval_f(coeffs_b, 3, tb, num_points_b, &waveform[num_points_a]); 

  KS_DESCRIPTION tmpdesc;
  ks_create_suffixed_description(tmpdesc, desc, "_cs%d", shift);

  status = ks_eval_wave(wave, tmpdesc, num_points_a + num_points_b, duration, waveform);
  KS_RAISE(status);

  if (wave->min_amp < 0) {
    return ks_error("%s: wave undershoots", __FUNCTION__);
  }

  ks_wave_multiplyval(wave, area / wave->area); /* normalize to account for round-off errors */
  ks_wave_compute_params(wave);

  if (wave->abs_max_amp > ks_syslimits_ampmax1(loggrd)) {
    return ks_error("%s: wave exceeds amplitude limit", __FUNCTION__);
  }
  if (wave->abs_max_slew > ks_syslimits_slewrate1(loggrd)) {
    return ks_error("%s: wave exceeds slew limit", __FUNCTION__);
  }


  return SUCCESS;

} /* ksdixon_eval_asym_spline_acq_wave() */

ksdixon_eval_trip_trap_acq_points()

STATUS ksdixon_eval_trip_trap_acq_points	(	float *	t,
		float *	G,
		int	shift,
		int	duration,
		float	area,
		float	Mc,
		float	g1,
		float	s
	)

Generates support points for an asymmetric piece-wise linear waveform with three bandwidths

[FSE Dixon with bandwidth-matched asymmetric readouts tailored for real-valued fat/water estimates, ISMRM 2023]

Parameters

[out]	t	6 times
[out]	G	6 amplitudes
[in]	shift	Desired CSE
[in]	duration	Duration (Will equal t[5] - t[0])
[in]	area	Required sample area
[in]	Mc	Mid plateau area
[in]	g1	Mid plateau amp
[in]	s	slewrate

Return values

STATUS

SUCCESS or FAILURE

                                            {

  float Ms = (area - Mc) / 2; /* Side area (area of 1st and 3rd trap) */

  /* Trap durations s.t. duration = T1 + T2 + T3 */
  float T2 = (Mc / g1);
  float center = duration/2 - abs(shift); /* Gradient echo (center of 2nd plateau) */
  float T1 = center - T2/2;
  float T3 = duration - T2 - T1;

  float t2_1   = - ( sqrt( s*T1*T1 + 2*g1*T1 - 2*Ms) - T1 * sqrt(s) ) / sqrt(s);
  float t2_2   =   ( sqrt( s*T1*T1 + 2*g1*T1 - 2*Ms) + T1 * sqrt(s) ) / sqrt(s);
  float t2 = t2_1 > 0 ? t2_1 : t2_2;

  float t1 = T1 - t2;     /* First trap plateau */
  float g0 = g1 + s * t2; /* First trap  amp */

  float t3_1   =   ( sqrt(s*T3*T3 - 2*T3*g1 + 2*Ms) + T3 * sqrt(s) ) / sqrt(s);
  float t3_2   = - ( sqrt(s*T3*T3 - 2*T3*g1 + 2*Ms) - T3 * sqrt(s) ) / sqrt(s);

  float t3 = t3_1 > T3 ? t3_2 : t3_1; /* Third trap ramp */
  float g2 = g1 - s*t3; /* Third trap amp */

  if (t3 < T3 &&
      t2 < T1 &&
      g0 > 0 &&
      g2 > 0) {
      if (shift <= 0) {
        t[0] = 0;             G[0] = g0;
        t[1] = t1;            G[1] = g0;
        t[2] = T1;            G[2] = g1;
        t[3] = T1+T2;         G[3] = g1;
        t[4] = T1+T2+t3;      G[4] = g2;
        t[5] = duration;      G[5] = g2;
      } else {
        t[5] = duration;                   G[5] = g0;
        t[4] = duration - (t1);            G[4] = g0;
        t[3] = duration - (T1);            G[3] = g1;
        t[2] = duration - (T1+T2);         G[2] = g1;
        t[1] = duration - (T1+T2+t3);      G[1] = g2;
        t[0] = duration - (duration);      G[0] = g2;
      }
      return SUCCESS;
  } else {
    return FAILURE;
  }

} /* ksdixon_eval_trip_trap_acq_points() */

ksdixon_dbw_ramptime()

float ksdixon_dbw_ramptime	(	const float	dur,
		const float	sample_area,
		const float	slewrate
	)

                                                                                           {
  const float root = sqrt(dur*dur - 4*sample_area/slewrate);
  if (dur < root) { /* Negative solution */
    ks_dbg("%s, negative solution - shouldn't happen.", __FUNCTION__);
    return ((dur + root)/2);
  } else {
    return ((dur - root)/2);
  }

} /* ksdixon_dbw_ramptime() */

ksdixon_eval_dbw_acq_waves()

STATUS ksdixon_eval_dbw_acq_waves	(	KS_WAVE *	lbw_wave,
		KS_WAVE *	hbw_wave,
		int *	lbw_cse,
		int *	hbw_cse,
		float *	pf,
		const float	min_pf,
		const int	duration,
		const float	area_without_pf,
		KS_DESCRIPTION	desc
	)

                                                                                                                                                                                                         {
  STATUS status;

  /* Brute-force optimize with free parameters nover (partial Fourier in read dir) and acq_duration of first trap */
  float max_nsa = -1.0f; /* target for optimization */

  int min_dur1 = RUP_FACTOR((int) (duration * min_pf), GRAD_UPDATE_TIME*2 );
  int dur_lbw;
  int dur_hbw;
  int opt_ramptime_lbw;
  float opt_ramptime_lbw_float;
  float cur_pf;
  const float slew = 0.97*ks_syslimits_slewrate1(loggrd);
  const float cse_at_acq_start = -duration / 2.0;
  const float dephasing_period = ksdixon_calc_dephasing_period(); /* usec */
  float best_time2echo_lbw = 0;
  float best_time2echo_hbw = 0;
  float cse[2] = {0.0, 0.0};

  float t_opt_lbw[2] = {0, 0};
  float G_opt_lbw[2] = {0, 0};
  float t_opt_hbw[4] = {0, 0, 0, 0};
  float G_opt_hbw[4] = {0, 0, 0, 0};

  int counter, opt_counter;
  for (dur_lbw = min_dur1; dur_lbw < duration; dur_lbw += GRAD_UPDATE_TIME ) {
    for (cur_pf = min_pf; cur_pf < 1; cur_pf += 0.01, counter++) {
      const float cur_sample_area = cur_pf * area_without_pf;
      /* Low bandwidth trapezoid. This excludes ramp up as that will be part of the crusher evaluation later */
      /* Ramptime if there was no raster */
      const float ramptime_lbw_float = ksdixon_dbw_ramptime(dur_lbw, cur_sample_area, slew);
      /* Snap ramptime to raster */
      int ramptime_lbw = RUP_FACTOR((int) ramptime_lbw_float, 2*GRAD_UPDATE_TIME); /* Why 2 times? because it's enforced in ks_eval_crusher */

      /* This shouldn't be necessary */
      ramptime_lbw += GRAD_UPDATE_TIME * 2;


      const float amp_lbw = cur_sample_area / (dur_lbw - ramptime_lbw);
      /* Since ramptime is now longer, it is necessary to adjust the slew rate to get the desired amp */
      const float slew_lbw_rastered = amp_lbw / ramptime_lbw;
      float t_lbw[4] = {      0,        0,  (float)(dur_lbw - ramptime_lbw),  (float)dur_lbw };
      float G_lbw[4] = {amp_lbw,  amp_lbw,                          amp_lbw,               0 };

      /* lbw is a late echo */
      float time2echo_lbw;
      status = ksdixon_dbw_time2echo_lbw(&time2echo_lbw, t_lbw, G_lbw, area_without_pf / 2.0);
      KS_RAISE(status);

      /* High bandwidth, use the remaining duration */

      /* 
                 rh
                  ____________
                 /|           |
                /<-----TOP--->|
               /__|___________|
              /|  |           |
             / |  |           |
            /  |  |           |
           /___|__|___________|
               rl
      */

      dur_hbw = duration - dur_lbw;
      /* Subtract the bottom rectangular area, making the ramptime solution in ksdixon_dbw_ramptime applicable */
      const float area_top = cur_sample_area - amp_lbw * (dur_hbw - ramptime_lbw);
      const float ramptime_hbw = ramptime_lbw + ksdixon_dbw_ramptime(dur_hbw, area_top, slew_lbw_rastered);
      if (isNaN(ramptime_hbw)) {
        continue;
      }
      const float amp_hbw = ramptime_hbw * slew_lbw_rastered;

      float t_hbw[4] = { 0, (float)ramptime_lbw, (float)ramptime_hbw, (float)dur_hbw };
      float G_hbw[4] = { 0,             amp_lbw,             amp_hbw,        amp_hbw };

      float time2echo_hbw;
      /* hbw is an early echo */
      status = ksdixon_dbw_time2echo_hbw(&time2echo_hbw, t_hbw, G_hbw, area_without_pf*(cur_pf - 0.5) + ramptime_lbw * amp_lbw / 2 );
      KS_RAISE(status);
      if (isNaN(time2echo_hbw)) {
        continue;
      }
      /* Gather chemical shift and evaluate the current gradient pair */
      cse[0] = cse_at_acq_start + time2echo_lbw;
      cse[1] = cse_at_acq_start + dur_lbw + time2echo_hbw;

      if (fabs(cse[0] - cse[1]) > dephasing_period) {
        /* If this failed, higher partial Fourier fractions will also fail */
        break;
      }

      float nsa = ksdixon_calc_nsa(cse[0], cse[1], 1.0 / (dur_lbw - ramptime_lbw), 1.0 / (dur_hbw - ramptime_lbw), dephasing_period);
      if (nsa > max_nsa) {
          max_nsa = nsa;
          opt_ramptime_lbw = ramptime_lbw;
          opt_ramptime_lbw_float = ramptime_lbw_float;
          t_opt_lbw[0] = 0;
          t_opt_lbw[1] = dur_lbw - ramptime_lbw;
          G_opt_lbw[0] = amp_lbw;
          G_opt_lbw[1] = amp_lbw;
          opt_counter = counter;
          t_opt_hbw[0] = 0;
          /* This point is necessary in order to amplitude match the two waves since ks_eval_coords2wave doesn't honor G[0] */
          t_opt_hbw[1] = GRAD_UPDATE_TIME;
          /* However, this alters the calculated times and area, so another correction is necessary */
          const float area_top_adj = area_top - GRAD_UPDATE_TIME * amp_lbw;
          const int rem_time = dur_hbw - ramptime_lbw - GRAD_UPDATE_TIME; /* On raster */
          const float ramptime_hbw_top = ksdixon_dbw_ramptime(rem_time, area_top_adj, slew_lbw_rastered);
          const float amp_top = area_top / (rem_time - ramptime_hbw_top / 2.0);
          t_opt_hbw[2] = GRAD_UPDATE_TIME + ramptime_hbw_top;
          t_opt_hbw[3] = dur_hbw - ramptime_lbw;
          G_opt_hbw[0] = amp_lbw;
          G_opt_hbw[1] = amp_lbw;
          G_opt_hbw[2] = amp_lbw + amp_top;
          G_opt_hbw[3] = amp_lbw + amp_top;
          best_time2echo_lbw = time2echo_lbw;
          best_time2echo_hbw = time2echo_hbw;
          *lbw_cse = (int)cse[0];
          *hbw_cse = (int)cse[1];
          *pf = cur_pf;
          status = SUCCESS;
      }
    } /* pf */
  } /* dur_lbw */
  ks_dbg("opt_counter = %d", opt_counter);
  if (status != SUCCESS || max_nsa < 0.0f) {
    *lbw_cse = 0;
    *hbw_cse = 0;
    return status;
  }

  KS_DESCRIPTION tmpdesc_lbw;
  ks_create_suffixed_description(tmpdesc_lbw, desc, "_lo_cs%d", *lbw_cse);
  status = ks_eval_coords2wave(lbw_wave, t_opt_lbw, G_opt_lbw, 2, GRAD_UPDATE_TIME, tmpdesc_lbw);
  KS_RAISE(status);
  KS_DESCRIPTION tmpdesc_hbw;
  ks_create_suffixed_description(tmpdesc_hbw, desc, "_hi_cs%d", *hbw_cse);
  status = ks_eval_coords2wave(hbw_wave, t_opt_hbw, G_opt_hbw, 4, GRAD_UPDATE_TIME, tmpdesc_hbw);
  KS_RAISE(status);


  ks_dbg("opt_ramptime_lbw = %d", opt_ramptime_lbw);
  ks_dbg("opt_ramptime_lbw_float = %.2f", opt_ramptime_lbw_float);
  ks_dbg("inner_ramparea = %.2f", opt_ramptime_lbw * G_opt_lbw[0] /2);
  ks_dbg("pf = %.2f", *pf);
  ks_dbg("area mismatch lbw = %.2f", lbw_wave->area - (*pf) * area_without_pf);
  ks_dbg("area mismatch hbw = %.2f", hbw_wave->area - (*pf) * area_without_pf);
  ks_dbg("sample_area = %.2f", (*pf) * area_without_pf);
  ks_dbg("lbw_wave_area = %.2f", lbw_wave->area);
  ks_dbg("hbw_wave_area = %.2f", hbw_wave->area);
  ks_dbg("best_time2echo_lbw = %.2f", best_time2echo_lbw);
  ks_dbg("best_time2echo_hbw = %.2f", best_time2echo_hbw);
  ks_dbg("cse = %d/%d", *lbw_cse, *hbw_cse);
  int i;
  for (i = 0; i < 2; i++) {
    ks_dbg("t_opt_lbw[%d] = %.2f", i, t_opt_lbw[i]);
    ks_dbg("G_opt_lbw[%d] = %.4f", i, G_opt_lbw[i]);
  }

  for (i = 0; i < 4; i++) {
    ks_dbg("t_opt_hbw[%d] = %.2f", i, t_opt_hbw[i]);
    ks_dbg("G_opt_hbw[%d] = %.4f", i, G_opt_hbw[i]);
  }
  return SUCCESS;

} /* ksdixon_eval_dbw_acq_waves() */

ksdixon_eval_flat_acq_wave()

STATUS ksdixon_eval_flat_acq_wave	(	KS_WAVE *	wave,
		int	duration,
		float	area,
		KS_DESCRIPTION	desc
	)

Generates an acq wave without chemical shift encoding (i.e. flat)

Parameters

[out]	wave	Acquisition wave
[in]	duration	Acq duration
[in]	area	Required sample area
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                {
  float t[2] = {0, (float)duration};
  float G[2] = {area/duration, area/duration};
  STATUS status = ks_eval_coords2wave(wave, t, G, 2, GRAD_UPDATE_TIME, desc);
  ks_create_suffixed_description(desc, desc, "_flat");
  return status;

} /* ksdixon_eval_flat_acq_wave() */

ksdixon_eval_inverted_ramp_acq_wave()

STATUS ksdixon_eval_inverted_ramp_acq_wave	(	KS_WAVE *	wave,
		int	shift,
		int	duration,
		float	area,
		KS_DESCRIPTION	desc
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

[out]	wave	ADDTEXTHERE
[in]	shift	ADDTEXTHERE
[in]	duration	ADDTEXTHERE
[in]	area	ADDTEXTHERE
[in]	desc	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                    {
  STATUS status;
  KS_DESCRIPTION tmpdesc;
  (void) tmpdesc;
  if (shift == 0) { /* Make a flat acq wave if in-phase echo is desired */
    return ksdixon_eval_flat_acq_wave(wave, duration, area, desc);
  }
  status = SUCCESS;
  return status;

} /* ksdixon_eval_inverted_ramp_acq_wave() */

ksdixon_eval_trip_trap_acq_wave()

STATUS ksdixon_eval_trip_trap_acq_wave	(	KS_WAVE *	wave,
		int	shift,
		int	duration,
		float	area,
		KS_DESCRIPTION	desc
	)

Wrapper to ksdixon_eval_trip_trap_acq_points() to get an acq wave (KS_WAVE) given a specific shift, duration, and area.

[FSE Dixon with bandwidth-matched asymmetric readouts tailored for real-valued fat/water estimates, ISMRM 2023]

Parameters

[out]	wave	Acquisition wave
[in]	shift	Desired CSE
[in]	duration	Acq duration
[in]	area	Required sample area
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                {
  STATUS status;
  KS_DESCRIPTION tmpdesc;
  if (shift == 0) { /* Make a flat acq wave if in-phase echo is desired */
    return ksdixon_eval_flat_acq_wave(wave, duration, area, desc);
  }

  if (duration % GRAD_UPDATE_TIME) {
    return ks_error("%s: duration (%d) must fit on the gradient raster (%d)", __FUNCTION__, duration, GRAD_UPDATE_TIME);
  }

  float t[6] = {KS_NOTSET};
  float G[6] = {KS_NOTSET};
  float amp = area / duration;
  float center_area = area * .25;
  status = ksdixon_eval_trip_trap_acq_points(t, G, shift, duration, area, center_area, amp, ks_syslimits_slewrate1(loggrd));
  KS_RAISE(status);
  ks_create_suffixed_description(tmpdesc, desc, "_cs%d", shift);
  status = ks_eval_coords2wave(wave, t, G, 6, GRAD_UPDATE_TIME, tmpdesc);

  return status;

} /* ksdixon_eval_trip_trap_acq_wave() */

ksdixon_amp_limit()

float ksdixon_amp_limit	(	float	board_amp_limit,
		float	fov
	)

Returns the highest amplitude allowed. This can either be the board max, or restricted by the FOV.

Parameters

[in]	board_amp_limit	[G/cm] Max board amplitude
[in]	fov	[mm] FOV

Returns

float

                                                          {
  float bw_amp_limit = ((float) (1.0 / (GAM * 2 * 1e-6) * (10.0 / fov)));
  float amp_limit = board_amp_limit < bw_amp_limit ? board_amp_limit : bw_amp_limit;

 return amp_limit;

} /* ksdixon_amp_limit() */

ksdixon_eval_inverted_ramps()

STATUS ksdixon_eval_inverted_ramps	(	KS_WAVE *	acq_waves,
		int *	shifts,
		int	points,
		const KS_KSPACE_DESIGN *	kdesign,
		KS_DESCRIPTION	desc
	)

                                                                                                                                      {
  STATUS status;
  KS_DESCRIPTION tmpdesc;
  /* get a resonable acq window and reference amplitude from the default trap */
  KS_READTRAP default_trap;
  status = ks_eval_design_readtrap(&default_trap, kdesign, 0, "default");
  KS_RAISE(status);

  float duration = default_trap.acq.duration;
  float sample_area = duration * default_trap.grad.amp;

  float slew = ks_syslimits_slewrate1(loggrd) * 0.95;
  float amp_limit = ksdixon_amp_limit(ks_syslimits_ampmax1(loggrd), kdesign->fov[XGRAD]);

  int i = 0;
  for (; i < points; i++) {
    sprintf(tmpdesc, "%s_invr_cs%d", desc, shifts[i]);
    if (shifts[i] == 0) {
      /* return a trap */
      float t[2] = {0, (float)duration};
      float G[2] = {default_trap.grad.amp, default_trap.grad.amp};
      status = ks_eval_coords2wave(&acq_waves[i], t, G, 2, GRAD_UPDATE_TIME, tmpdesc);
      KS_RAISE(status);
    } else {
      float t[4] = {0};
      float G[4] = {0};
      status = ksdixon_eval_inverted_ramp_acq_points(t, G, shifts[i], duration, sample_area, slew, amp_limit);
      KS_RAISE(status);
      status = ks_eval_coords2wave(&acq_waves[i], t, G, 4, GRAD_UPDATE_TIME, tmpdesc);
      KS_RAISE(status);
    }
  }
  return SUCCESS;

} /* ksdixon_eval_inverted_ramps() */

ksdixon_eval_inverted_ramp_acq_points()

STATUS ksdixon_eval_inverted_ramp_acq_points	(	float *	t,
		float *	G,
		int	shift,
		int	duration,
		float	area,
		float	s,
		float	amp_limit
	)

Generates support points for an inverted ramp waveform.

Parameters

[out]	t	4 times
[out]	G	4 amplitudes
[in]	shift	Desired CSE
[in]	duration	Duration (Will equal t[3] - t[0])
[in]	area	Required sample area
[in]	s	slewrate
[in]	amp_limit	Maximum amplitude

                                                      {

  float plat_amp = area / (2 * fabs(shift) + duration);
  float peak_amp = plat_amp + 2 * sqrt(plat_amp * fabs(shift) * s);
  float ramp_amp = (peak_amp - plat_amp);

  if (peak_amp > amp_limit) {
    peak_amp = amp_limit;
    ramp_amp = (peak_amp - plat_amp);

    float g0 = peak_amp;
    float g1 = area / (duration + 2*fabs(shift) );
    float T2 = (g0-g1) / s;
    float T1 = (area - T2*(g0 + g1)/2 - (duration - T2) * g1) / (g0-g1);
    if (T1 < 0) {
      return ks_error("%s: Negative initial plateau. Something is wrong.", __FUNCTION__);
    }
    float T3 = duration - T2 - T1;
    if (shift < 0) {
      G[0] = g0;
      G[1] = g0;
      G[2] = g1;
      G[3] = g1;
      t[0] = 0;
      t[1] = T1;
      t[2] = T1 + T2;
      t[3] = duration;
    } else {
      G[0] = g1;
      G[1] = g1;
      G[2] = g0;
      G[3] = g0;
      t[0] = 0;
      t[1] = T3;
      t[2] = T3 + T2;
      t[3] = duration;
    }

    return SUCCESS;
  }

  if ((int)(ramp_amp / s) > (int)(duration / 2 - fabs(shift))) {
    return ks_error("Too large a shift for the slew limit!");
  }

  if (shift > 0) {
    G[0] = plat_amp;
    t[0] = 0.0;
    G[1] = plat_amp;
    t[1] = duration - ramp_amp / s;
    G[2] = peak_amp;
    t[2] = duration;
  } else if (shift == 0) {
    G[0] = plat_amp;
    t[0] = 0.0;
    t[2] = plat_amp;
    t[2] = duration;
    G[1] = plat_amp;
    t[1] = duration/2;
  } else { /* negative shift */
    G[0] = peak_amp;
    t[0] = 0.0;
    G[1] = plat_amp;
    t[1] = ramp_amp / s;
    G[2] = plat_amp;
    t[2] = duration; 
  }
  t[3] = t[2];
  G[3] = G[2];

  return SUCCESS;

} /* ksdixon_eval_inverted_ramp_acq_points() */

ksdixon_init_settings()

void ksdixon_init_settings

(

KSDIXON_DESIGN *const

dixon

)

Resets a KSDIXON_DESIGN to KSDIXON_INIT_DESIGN

Parameters

dixon

The KSDIXON_DESIGN to reset

                                                         {
  KSDIXON_DESIGN default_settings = KSDIXON_INIT_DESIGN;

  *dixon = default_settings;

} /* ksdixon_init_settings() */

ksdixon_setup()

STATUS ksdixon_setup	(	KSDIXON_DESIGN *	pdixon,
		const int	opnum_dixon_mode,
		_cvfloat *	_opuser_dixon_mode,
		const int	opnum_dixon_tr_mode,
		_cvfloat *	_opuser_dixon_tr_mode,
		const int	opnum_dixon_shift,
		_cvfloat *	_opuser_dixon_shift
	)

ksdixon_setup

Intended to be called using the macro KSDIXON_SETUP:

Implementation e file

  // Select opuser numbers for dixon options
  @global
  #define OPUSER_NUM_DIXON_MODE 17
  #define OPUSER_NUM_DIXON_TR_MODE 18
  #define OPUSER_NUM_DIXON_SHIFT 19

ksfse2_setup_stretched():
  // Fill ksfse2_dixon_settings with dixon related values.
  status = KSDIXON_SETUP(&ksfse2_dixon_settings, OPUSER_NUM_DIXON_MODE, OPUSER_NUM_DIXON_TR_MODE, OPUSER_NUM_DIXON_SHIFT);

                                                                                 {
  (void)opnum_dixon_mode;
  (void)opnum_dixon_tr_mode;
  (void)opnum_dixon_shift;

  STATUS status;
  ksdixon_init_settings(pdixon);

  pdixon->mode = (KSDIXON_MODE) *_opuser_dixon_mode->addr;

  pdixon->points = 2;

  if (pdixon->mode != KSDIXON_DBW) {

    pdixon->tr_mode = (KSDIXON_TR_MODE) *_opuser_dixon_tr_mode->addr;
    pdixon->shifts[0] = 0;
    pdixon->shifts[1] = RUP_GRD((int) *_opuser_dixon_shift->addr);

  }

  status = ksdixon_eval_validate_settings(pdixon);
  KS_RAISE(status);

  return SUCCESS;

} /* ksdixon_setup() */

ksdixon_eval_validate_settings()

STATUS ksdixon_eval_validate_settings

(

const KSDIXON_DESIGN *

dixon

)

Validates KSDIXON_DESIGN to make sure settings are compatible. Right now it only checks if DBW is in single TR mode.

Parameters

[in]

dixon

The settings

Return values

STATUS

SUCCESS or FAILURE

                                                                   {
  if ((dixon->tr_mode == KSDIXON_MULTI_TR) && (dixon->mode == KSDIXON_DBW)) {
    return ks_error("%s: Multi-TR not supported for dBW readouts", __FUNCTION__);
  }
  return SUCCESS;

} /* ksdixon_eval_validate_settings() */

ksdixon_update_UI()

void ksdixon_update_UI	(	const KSDIXON_MODE	dixon_mode,
		const KS_KSPACE_DESIGN *const	kdesign,
		const int	opnum_dixon_mode,
		_cvfloat *	_opuser_dixon_mode,
		const int	opnum_dixon_tr_mode,
		_cvfloat *	_opuser_dixon_tr_mode,
		const int	opnum_dixon_shift,
		_cvfloat *	_opuser_dixon_shift
	)

ksdixon_update_UI

Intended to be called using the macro KSDIXON_UPDATE_UI:

Currently not used. The idea is to update dixon opusers according to kdesign and dixon_mode.

                                                                                   {
  (void)_opuser_dixon_mode;
  (void)_opuser_dixon_tr_mode;
  (void)_opuser_dixon_shift;
  (void)kdesign;

  switch (dixon_mode) {
    case KSDIXON_OFF: {
      piuset |= (1 << opnum_dixon_mode);
      piuset &= ~( (1 << opnum_dixon_tr_mode) | (1 << opnum_dixon_shift) );
      break;
    }
    case KSDIXON_SHIFTED:
    case KSDIXON_ASYM_SPLINE:
    case KSDIXON_ASYM_INVRAMP:
    case KSDIXON_TRIP_TRAP: {
      piuset |= (1 << opnum_dixon_tr_mode);
      piuset |= (1 << opnum_dixon_shift);
      break;
    }
    case KSDIXON_DBW:
    case KSDIXON_DBW_STAGGERED: {
      piuset &= ~( (1 << opnum_dixon_tr_mode) | (1 << opnum_dixon_shift) );
      break;
    }
    default:
      break;
  }

} /* ksdixon_update_UI() */

ksdixon_init_UI()

void ksdixon_init_UI	(	float	rbw,
		int	xres,
		float	fov,
		const int	opnum_dixon_mode,
		_cvfloat *	_opuser_dixon_mode,
		const int	opnum_dixon_tr_mode,
		_cvfloat *	_opuser_dixon_tr_mode,
		const int	opnum_dixon_shift,
		_cvfloat *	_opuser_dixon_shift
	)

See documentation for macro KSDIXON_INIT_UI.

                                                                                 {

  piuset |= (1 << opnum_dixon_mode);
  _opuser_dixon_mode->minval = 0;
  _opuser_dixon_mode->maxval = KSDIXON_MODE_NUM_MODES-1;
  _opuser_dixon_mode->defval = KSDIXON_OFF;
  _opuser_dixon_mode->descr = KSDIXON_MODE_DESCR;
  KSDIXON_MODE mode = (KSDIXON_MODE) *_opuser_dixon_mode->addr;

  switch (mode) {
    case KSDIXON_OFF: {
      piuset &= ~( (1 << opnum_dixon_tr_mode) | (1 << opnum_dixon_shift) );
      break;
    }
    case KSDIXON_SHIFTED:
    case KSDIXON_ASYM_SPLINE:
    case KSDIXON_ASYM_INVRAMP:
    case KSDIXON_TRIP_TRAP: {
      piuset |= (1 << opnum_dixon_tr_mode);
      _opuser_dixon_tr_mode->minval = 0;
      _opuser_dixon_tr_mode->maxval = 1;
      _opuser_dixon_tr_mode->defval = KSDIXON_SINGLE_TR;
      _opuser_dixon_tr_mode->descr = KSDIXON_TR_MODE_DESCR;

      piuset |= (1 << opnum_dixon_shift);
      int tsp = ks_calc_bw2tsp(rbw);
      int max_shift = ksdixon_max_shift(mode, xres * tsp, xres, fov);
      _opuser_dixon_shift->minval = -5000;
      _opuser_dixon_shift->maxval = 5000;
      _opuser_dixon_shift->defval = 0;
      KS_DESCRIPTION tmpdesc;
      sprintf(tmpdesc,"%s (max %d)",KSDIXON_SHIFT_DESC, max_shift);
      _opuser_dixon_shift->descr = strdup(tmpdesc);
      break;
    }
    case KSDIXON_DBW:
    case KSDIXON_DBW_STAGGERED: {
      piuset &= ~( (1 << opnum_dixon_tr_mode) | (1 << opnum_dixon_shift) );
      break;
    }
    default:
      piuset &= ~(1 << opnum_dixon_tr_mode);
      piuset &= ~(1 << opnum_dixon_shift);
      break;
  }

} /* ksdixon_init_UI() */

ksdixon_eval_setup_gre_readouts()

STATUS ksdixon_eval_setup_gre_readouts	(	KS_ECHOTRAIN *	echotrain,
		KSDIXON_DESIGN *	dixon,
		const KS_KSPACE_DESIGN *	kdesign,
		const int	acq_start_time,
		const int	TR_index,
		KS_DESCRIPTION	desc
	)

Sets up GRE readouts. Currently only supports KSDIXON_DBW

Parameters

[out]	echotrain	Pointer to a KS_ECHOTRAIN, which keeps track of states and what to pulsegen
[in,out]	dixon	Dixon settings
[in]	kdesign	K-space deisng
[in]	acq_start_time	Time when first acquisition is turned on
[in]	TR_index	TR index for multi-TR support (currently not used as only DBW is supported)
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                                                           {
  STATUS status;
  int readout;

  if (dixon->mode == KSDIXON_OFF) {
    return ks_error("%s should not be called for mode KSDIXON_OFF", __FUNCTION__);
  }

  switch (dixon->mode) { 

    case KSDIXON_DBW:
      {
        if (TR_index > 0) {
          return ks_error("%s: Multi-TR not supported for dual BW Dixon", __FUNCTION__);
        }
        if (dixon->points != 2) {
          return ks_error("%s: dual BW Dixon requires exactly two points, not %d", __FUNCTION__, dixon->points);
        }
        KS_READWAVE* trap1 = &echotrain->readwaves[0];
        KS_READWAVE* trap2 = &echotrain->readwaves[1];
        status = ksdixon_eval_gre_dualBW_readwave(trap1, trap2, kdesign, dixon, acq_start_time, desc);
        KS_RAISE(status);
        for(readout = 0; readout < 2; readout++) {
          echotrain->controls[readout].state[0].kspace_index = readout;
          echotrain->controls[readout].pg.object_index = readout;
          echotrain->controls[readout].pg.read_type = KS_READ_WAVE;
          echotrain->controls[readout].pg.ampscale = (readout==0) ? 1.0f : -1.0f;
        }
        return SUCCESS;
      }

    default:
      return ks_error("%s: KSDIXON_MODE %d not implemented", __FUNCTION__, dixon->mode);
  }

  return SUCCESS;

} /* ksdixon_eval_setup_gre_readouts() */

ksdixon_set_design()

STATUS ksdixon_set_design	(	KSDIXON_DESIGN *const	dixon_design,
		const KS_KSPACE_DESIGN *const	kdesign
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

dixon_design	ADDTEXTHERE
kdesign	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                     {
  STATUS status;
  /* get a resonable acq window from the default trap */
  KS_READTRAP default_trap;
  KS_KSPACE_DESIGN default_kdesign = *kdesign;
  status = ks_eval_design_readtrap(&default_trap, &default_kdesign, 0.0, "default");
  KS_RAISE(status);
  dixon_design->acq_duration = RUP_FACTOR(default_trap.acq.duration, 2*GRAD_UPDATE_TIME);
  dixon_design->sample_area = dixon_design->acq_duration * default_trap.grad.amp;
  dixon_design->area2center = kdesign->nover[XGRAD] == 0 ? dixon_design->sample_area / 2
                                                         : dixon_design->sample_area * 1.0f/(1 + kdesign->res[XGRAD] / (2.0f * kdesign->nover[XGRAD]) ); /* Assume early echo */

  return SUCCESS;

} /* ksdixon_set_design() */

ksdixon_setup_fse_echotrain()

STATUS ksdixon_setup_fse_echotrain	(	KS_ECHOTRAIN_DESIGN *const	echotrain_design,
		KSDIXON_STATE *const	dixon_state,
		const KSDIXON_DESIGN *const	dixon_design,
		const int	ETL,
		float	xcrusher_area,
		const KS_DESCRIPTION	description
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	echotrain_design	ADDTEXTHERE
	dixon_state	ADDTEXTHERE
	dixon_design	ADDTEXTHERE
[in]	ETL	ADDTEXTHERE
[in]	xcrusher_area	ADDTEXTHERE
[in]	description	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                     {
ks_enum_crusher_strategy crusher_strategy;
STATUS status;
int i;
int cse;
int spin_echo;

switch (dixon_design->mode) {
  case KSDIXON_OFF:
    return SUCCESS;
  case KSDIXON_SHIFTED:
    return SUCCESS;
  case KSDIXON_DBW_STAGGERED:
  case KSDIXON_DBW: {
    KS_DESCRIPTION tmpdesc;
    ks_create_suffixed_description(tmpdesc, description, "_dbw");

    /* Allocate some memory for the acq waves */
    ks_init_design_readwave(&echotrain_design->readwave_designs[0]);
    ks_init_design_readwave(&echotrain_design->readwave_designs[1]);
    int numstates_per_readwaves = dixon_design->mode == KSDIXON_DBW_STAGGERED ? 2 : 1;
    for (i = 0; i < 2; i++) {
      KS_WAVE* acq_ptr = (KS_WAVE*) realloc(echotrain_design->readwave_designs[i].acq_waves, sizeof(KS_WAVE) * numstates_per_readwaves);
      if (!acq_ptr) {
        return KS_THROW("Reallocation of dbw acq wave #%d failed", i);
      }
      echotrain_design->readwave_designs[i].acq_waves = acq_ptr;
    }

    float min_pfx = (dixon_design->min_pfx >= 0.5) ? dixon_design->min_pfx
                                                   : 1;
    float pfx = 0;

    /* Generate dbw waves */
    KS_WAVE* lbw_waves = echotrain_design->readwave_designs[0].acq_waves;
    KS_WAVE* hbw_waves = echotrain_design->readwave_designs[1].acq_waves;
    ksdixon_eval_dbw_acq_waves(lbw_waves,
                               hbw_waves,
                               &dixon_state->shifts[0],
                               &dixon_state->shifts[1],
                               &pfx,
                               min_pfx,
                               dixon_design->acq_duration,
                               dixon_design->sample_area,
                               tmpdesc);
    int inner_ramp = (dixon_design->acq_duration - lbw_waves->duration - hbw_waves->duration) / 2;
    if (inner_ramp % GRAD_UPDATE_TIME) {
      ks_error("%s: dbw inner ramps (duration %d) are not gradient raster compatible", __FUNCTION__, inner_ramp);
      KS_RAISE(FAILURE);
    }
    /* Waves are centered around spin echo so pg_shifts must correct for that */
    int pg_shifts[2] = {(lbw_waves->duration - dixon_design->acq_duration)/2,
                        (hbw_waves->duration - dixon_design->acq_duration)/2 + lbw_waves->duration + 2*inner_ramp};
    dixon_state->deadtime = 2 * inner_ramp;
    ks_dbg("DBW inner ramps are %d", inner_ramp);

    /* Setup readwave designs */
    echotrain_design->readwave_designs[0].flag_symmetric_padding = 0;
    echotrain_design->readwave_designs[0].pf = -pfx; /* Negative because it is a late echo */
    echotrain_design->readwave_designs[0].numstates = numstates_per_readwaves;
    echotrain_design->readwave_designs[0].rbw = ks_calc_tsp2bw(2);

    echotrain_design->readwave_designs[0].pre_crusher.area = xcrusher_area;
    echotrain_design->readwave_designs[0].pre_crusher.strategy = KS_CRUSHER_STRATEGY_PLATEAU;
    echotrain_design->readwave_designs[0].pre_crusher.slew = 0; /* No preference */
    echotrain_design->readwave_designs[0].post_crusher.area = 0; /* Inner ramps */
    echotrain_design->readwave_designs[0].post_crusher.strategy = KS_CRUSHER_STRATEGY_FASTEST;
    echotrain_design->readwave_designs[0].post_crusher.slew = lbw_waves->abs_max_slew; /* Inner ramps must match */


    echotrain_design->readwave_designs[1].flag_symmetric_padding = 0;
    echotrain_design->readwave_designs[1].pf = pfx;
    echotrain_design->readwave_designs[1].numstates = numstates_per_readwaves; 
    echotrain_design->readwave_designs[1].rbw = ks_calc_tsp2bw(2);

    echotrain_design->readwave_designs[1].pre_crusher.area = 0; /* Inner ramp  */
    echotrain_design->readwave_designs[1].pre_crusher.strategy = KS_CRUSHER_STRATEGY_FASTEST;
    echotrain_design->readwave_designs[1].pre_crusher.slew = echotrain_design->readwave_designs[0].post_crusher.slew; /* Inner ramps must match */
    echotrain_design->readwave_designs[1].post_crusher.area = xcrusher_area;
    echotrain_design->readwave_designs[1].post_crusher.strategy = KS_CRUSHER_STRATEGY_PLATEAU;
    echotrain_design->readwave_designs[1].post_crusher.slew = 0; /* No preference */


    /*negated state of each acq wave */
    /* lbw */
    lbw_waves[1] = lbw_waves[0];
    ks_waveform_multiplyval(lbw_waves[1].waveform,  -1, lbw_waves[1].res);

    /* hbw */
    hbw_waves[1] = hbw_waves[0];
    ks_waveform_multiplyval(hbw_waves[0].waveform,  -1, hbw_waves[0].res);
    ks_wave_compute_params(&hbw_waves[0]);

    /* Setup read controls */
    for(spin_echo=0; spin_echo < ETL; spin_echo++) {
      for(cse=0; cse < 2; cse++) {
          i = spin_echo * 2 + cse;
          echotrain_design->controls[i].state[0].kspace_index = cse;
          echotrain_design->controls[i].state[0].wave_state = 0;
          echotrain_design->controls[i].pg.object_index = cse;
          echotrain_design->controls[i].pg.read_type = KS_READ_WAVE;
          echotrain_design->controls[i].pg.shift = pg_shifts[cse];
          echotrain_design->controls[i].pg.ampscale = 1.0f; /* PG both waves positive and change their state in scan */
      }
    }
    break;
  }
  case KSDIXON_TRIP_TRAP:
  case KSDIXON_ASYM_INVRAMP:
  case KSDIXON_ASYM_SPLINE: {
    if (dixon_design->points != 2) {
      return ks_error("%s: KSDIXON_ASYM_SPLINE only implemented for two points (%d were prescribed)", __FUNCTION__, dixon_design->points);
    }

    ks_init_design_readwave(&echotrain_design->readwave_designs[0]);

    /* Generate the acq waves */
    for (i = 0; i < dixon_design->points; i++) {
      ks_init_design_readwave(&echotrain_design->readwave_designs[i]);
      KS_WAVE* acq_ptr = NULL;
      if (dixon_design->tr_mode == KSDIXON_MULTI_TR) {
        if (i == 0) {
          acq_ptr = (KS_WAVE*) realloc(echotrain_design->readwave_designs[0].acq_waves, dixon_design->points * sizeof(KS_WAVE));
          if (!acq_ptr) {
            return KS_THROW("Reallocation of %d acq waves failed", dixon_design->points);
          }
          echotrain_design->readwave_designs[0].acq_waves = acq_ptr;
        } else {
          acq_ptr = &echotrain_design->readwave_designs[0].acq_waves[i];
        }
      } else if (dixon_design->tr_mode == KSDIXON_SINGLE_TR) {
        acq_ptr = (KS_WAVE*) realloc(echotrain_design->readwave_designs[i].acq_waves, sizeof(KS_WAVE));
        if (!acq_ptr) {
          return KS_THROW("Reallocation of acq wave failed");
        }
        echotrain_design->readwave_designs[i].acq_waves = acq_ptr;
      }


      KS_DESCRIPTION tmpdesc;

      if (dixon_design->mode == KSDIXON_ASYM_SPLINE) {
        ks_create_suffixed_description(tmpdesc, description, "_spline%d", i);
        status = ksdixon_eval_asym_spline_acq_wave(acq_ptr, dixon_design->shifts[i], dixon_design->acq_duration, dixon_design->sample_area, tmpdesc);
        dixon_state->shifts[i] = dixon_design->shifts[i];
        crusher_strategy = KS_CRUSHER_STRATEGY_PLATEAU;
      } else if (dixon_design->mode == KSDIXON_TRIP_TRAP) {
        ks_create_suffixed_description(tmpdesc, description, "_pwl%d", i);
        status = ksdixon_eval_trip_trap_acq_wave(acq_ptr, dixon_design->shifts[i], dixon_design->acq_duration, dixon_design->sample_area, tmpdesc);
        dixon_state->shifts[i] = dixon_design->shifts[i];
        crusher_strategy = KS_CRUSHER_STRATEGY_FASTEST;
      } else if (dixon_design->mode == KSDIXON_ASYM_INVRAMP) {
        ks_create_suffixed_description(tmpdesc, description, "_invr%d", i);
        status = ksdixon_eval_inverted_ramp_acq_wave(acq_ptr, dixon_design->shifts[i], dixon_design->acq_duration, dixon_design->sample_area, tmpdesc);
        dixon_state->shifts[i] = dixon_design->shifts[i];
        crusher_strategy = KS_CRUSHER_STRATEGY_FASTEST;
      }
      KS_RAISE(status);
    }

    /* Setup the readwave designs and readout controls */
    echotrain_design->readwave_designs[0].flag_symmetric_padding = 1;
    echotrain_design->readwave_designs[0].pf = 1.0; /* No pf */
    echotrain_design->readwave_designs[0].pre_crusher.area = xcrusher_area; /* Changed later for single TR */
    echotrain_design->readwave_designs[0].post_crusher.area = xcrusher_area; /* Changed later for single TR */
    echotrain_design->readwave_designs[0].pre_crusher.strategy = crusher_strategy;
    echotrain_design->readwave_designs[0].post_crusher.strategy = crusher_strategy;
    echotrain_design->readwave_designs[0].pre_crusher.slew = 0; /* No preference */
    echotrain_design->readwave_designs[0].post_crusher.slew = 0; /* No preference */

    /* TODO: Throw if not set */
    echotrain_design->readwave_designs[0].pre_crusher.ampmax = dixon_design->amp_crusher;
    echotrain_design->readwave_designs[0].post_crusher.ampmax = dixon_design->amp_crusher;
    echotrain_design->readwave_designs[0].pre_crusher.slew = dixon_design->slew_crusher;
    echotrain_design->readwave_designs[0].post_crusher.slew = dixon_design->slew_crusher;

    if (dixon_design->tr_mode == KSDIXON_SINGLE_TR) {
      /* In single TR mode, each readwave is stored separately with a single state */
      for (i = 0; i < dixon_design->points; i++) {
        KS_WAVE* acq_ptr = echotrain_design->readwave_designs[i].acq_waves;
        echotrain_design->readwave_designs[i] = echotrain_design->readwave_designs[0];
        echotrain_design->readwave_designs[i].acq_waves = acq_ptr;
        echotrain_design->readwave_designs[i].numstates = 1; /* States per readwave_design */
      }

      /* It's important that crusher area matches or we get stimulated echoesin FSE. We figure out the max area out first */
      KS_WAVE tmp_crusher = KS_INIT_WAVE;
      ks_crusher_constraints crusher_constraints = echotrain_design->readwave_designs[0].pre_crusher;
      float max_crusher_area = xcrusher_area;
      for (i = 0; i < dixon_design->points; i++) {
        KS_WAVE* acq_ptr = echotrain_design->readwave_designs[i].acq_waves;

        /* Estimate coords */
        const float start_amp = 1.5f * acq_ptr->waveform[0] - 0.5f * acq_ptr->waveform[1];
        const float end_amp = 1.5f * acq_ptr->waveform[acq_ptr->res - 1] - 0.5f * acq_ptr->waveform[acq_ptr->res - 2];

        ks_eval_crusher(&tmp_crusher, start_amp, crusher_constraints);
        if (tmp_crusher.area > max_crusher_area) {
          max_crusher_area = tmp_crusher.area;
        }

        ks_eval_crusher(&tmp_crusher, end_amp, crusher_constraints);
        if (tmp_crusher.area > max_crusher_area) {
          max_crusher_area = tmp_crusher.area;
        }
      }

      echotrain_design->readwave_designs[0].pre_crusher.area =  max_crusher_area;
      echotrain_design->readwave_designs[0].post_crusher.area = max_crusher_area;


      for(spin_echo=0; spin_echo < ETL; spin_echo++) {
        cse = spin_echo % dixon_design->points;
        echotrain_design->controls[spin_echo].state[0].kspace_index = cse; /* CSEs should be stored in separate k-spaces */
        echotrain_design->controls[spin_echo].state[0].wave_state = 0;
        echotrain_design->controls[spin_echo].pg.object_index = cse; /* So many readwaves */
        echotrain_design->controls[spin_echo].pg.read_type = KS_READ_WAVE;
        echotrain_design->controls[spin_echo].pg.ampscale = 1;
      }
    } else if (dixon_design->tr_mode == KSDIXON_MULTI_TR) {
      /* There is only one readwave with multiple states in multi TR mode */
      echotrain_design->readwave_designs[0].numstates = dixon_design->points;

      /* map cse to wavestate and kspace_index */
      for(spin_echo=0; spin_echo < ETL; spin_echo++) { /* readout loop */
        echotrain_design->controls[spin_echo].pg.read_type = KS_READ_WAVE;
        echotrain_design->controls[spin_echo].pg.ampscale = 1;
        echotrain_design->controls[spin_echo].pg.object_index = 0; /* As mentioned above, there is only one readwave in this case */
        for (cse = 0; cse < dixon_design->points; cse++) { /* state loop */
          echotrain_design->controls[spin_echo].state[cse].wave_state = cse;
          echotrain_design->controls[spin_echo].state[cse].kspace_index = cse;
        }
      }
    }
    break;
  }
  default:
    return KS_THROW("Dixon mode %d has not yet been implemented with the new latest and greatest rewrite.", dixon_design->mode);
    break;
  }

  return SUCCESS;
}

ksdixon_eval_trip_traps()

STATUS ksdixon_eval_trip_traps	(	KS_WAVE *	acq_waves,
		int *	shifts,
		int	points,
		const KS_KSPACE_DESIGN *	kdesign,
		KS_DESCRIPTION	desc
	)

Generates points readwaves, each has PWL3 acq waveform and trapezoidal crushers

Parameters

[out]	acq_waves	Output readwaves. Array with points elements must be allocated before
[in]	shifts	[us] Array of CSE
[in]	points	Number of readwaves to generate
[in]	kdesign	A KS_KSPACE_DESIGN, required to calculate sample area and acq duration
[in]	desc	Description

Return values

STATUS

SUCCESS or FAILURE

Variable Documentation

piuset

int piuset