ksepg.cc File Reference

#include <ksepg.h>
#include <stdlib.h>
#include <math.h>

int	rhkacq_uid
int	ks_plot_filefmt
int	ks_plot_kstmp
char	ks_psdname [256]
KSEPG_RELAXATION	ksepg_create_relaxation (double time_delta, double T1, double T2)
KSEPG_STATE	ksepg_create_state (int size)
void	ksepg_destruct_state (KSEPG_STATE *state)
KSEPG_STATE	ksepg_duplicate_state (const KSEPG_STATE *state)
void	ksepg_dephase (KSEPG_STATE *state, KSEPG_RELAXATION rel_op)
void	ksepg_compute_rotation (KS_MAT3x3 *rot, double flip)
void	ksepg_apply_rotation (KSEPG_STATE *state, KS_MAT3x3 *rot_p)
void	ksepg_refocus (KSEPG_STATE *state, KS_MAT3x3 *rot)
void	ksepg_next_echo (KSEPG_STATE *state, KS_MAT3x3 *rot, KSEPG_RELAXATION rel_op)
STATUS	ksepg_get_intensities (const double *flip_angles, int num_echoes, KSEPG_STATE state, double *signal_intensities, double echo_spacing, double T1, double T2)
STATUS	ksepg_get_intensities_std (const double *flip_angles, double *signal_intensities, int num_echoes, double echo_spacing, double T1, double T2)
STATUS	ksepg_get_intensities_long (const double *flip_angles, double *signal_intensities, int num_echoes, double echo_spacing, double T1, double T2)
double	ksepg_solve_trig_eq (double a, double b, double c)
double	ksepg_track_step (double Fp, double Fm, double Z, double target_Fm, KSEPG_RELAXATION rel_op)
STATUS	ksepg_track_intensity_curve (double *flip_angles, int num_echoes, KSEPG_STATE state, double scaling_factor, const double *target_curve, double echo_spacing, double T1, double T2)
STATUS	ksepg_track_intensity_std (double *flip_angles, int num_echoes, double target_intensity, double echo_spacing, double T1, double T2)
STATUS	ksepg_track_intensity_long (double *flip_angles, int num_echoes, double target_intensity, double echo_spacing, double T1, double T2)
void	ks_print_epg_sim (const float *flip_angles, const double *si_rel, const double *si_norel, const double *effTEs, const int num_echoes, const double echo_spacing, const double T1, const double T2, KS_DESCRIPTION description)
void	ks_predownload_plot_epg_sim (KS_DESCRIPTION description)
STATUS	ksepg_calc_effTE (const float *flip_angles, double *effTEs, int num_echoes, double echo_spacing, double T1, double T2, KS_DESCRIPTION description)
double	ksepg_max_angle (double *flip_angles, int num_echoes)
STATUS	ksepg_optimize_constant (double *intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_constant_std (double *intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2)
STATUS	ksepg_optimize_constant_long (double *intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2)
STATUS	ksepg_generate_intensity_triangle_impl (double *target_intensities, int num_echoes, double center_position, double min_intensity, double max_intensity, double half_size)
STATUS	ksepg_generate_intensity_periodic_triangle (double *target_intensities, int num_echoes, double center_position, double min_intensity, double max_intensity)
STATUS	ksepg_generate_intensity_cropped_triangle (double *target_intensities, int num_echoes, double center_position, double min_intensity, double max_intensity)
STATUS	ksepg_generate_intensity_ramp (double *target_intensities, int num_echoes, double center_position, double min_intensity, double max_intensity)
STATUS	ksepg_generate_intensity_linear (double *target_intensities, int num_echoes, double center_position, double min_intensity, double max_intensity)
STATUS	ksepg_track_bounded_curve (double *flip_angles, double *target_intensities, int num_echoes, double min_intensity, double max_intensity, double echo_spacing, double T1, double T2, double center_position, STATUS(*generate_bounded_curve)(double *, int, double, double, double), const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_bounded_curve (double *min_intensity, double *max_intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, double center_position, STATUS(*generate_bounded_curve)(double *, int, double, double, double), const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_bounded_curve_std (double *min_intensity, double *max_intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, double center_position, STATUS(*generate_bounded_curve)(double *, int, double, double, double))
STATUS	ksepg_optimize_bounded_curve_long (double *min_intensity, double *max_intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, double center_position, STATUS(*generate_bounded_curve)(double *, int, double, double, double))
STATUS	ksepg_optimize_rampdown (double *end_intensity, double start_intensity, double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, const KSEPG_STATE starting_state)
STATUS	ksepg_half_intensity (double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, KSEPG_TAIL tail_type)
STATUS	ksepg_optimize_expdown (double *flip_angles, int num_echoes, double max_flip, double echo_spacing, double T1, double T2, const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_half_exponential (double *flip_angles, int num_echoes, double max_flip)
STATUS	ksepg_optimize_half_free (double *flip_angles, int num_echoes, double max_flip)
STATUS	ksepg_optimize_linear_free (double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2, const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_linear_free_std (double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2)
STATUS	ksepg_optimize_linear_free_long (double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2)
STATUS	ksepg_optimize_linear_free_center_echo (int *_center_echo, int target_eff_TE, int linear_sweep_center_echo, double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2, const KSEPG_STATE starting_state)
STATUS	ksepg_optimize_linear_free_center_echo_std (int *center_echo, int target_eff_TE, int linear_sweep_center_echo, double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2)
STATUS	ksepg_optimize_linear_free_center_echo_long (int *center_echo, int target_eff_TE, int linear_sweep_center_echo, double *flip_angles, int num_echoes, double max_flip, double min_intensity_weight, double echo_spacing, double T1, double T2)
STATUS	ksepg_sineTRAPS (double *flip_angles, int etl, int n23, int n4, double an1, double an23, double an4)
STATUS	ksepg_sineTRAPS2 (float *flip_angles, int etl, KSEPG_TRAPS_DESIGN traps_design, float scaling_factor)
double	ksepg_calc_relSAR (double *flip_angles, int etl)

Detailed Description

This file contains sequence-independent implementations for extended phase graphs

Feature module: Extended Phase Graphs (ksepg.[h,cc])

Function Documentation

ksepg_create_relaxation()

KSEPG_RELAXATION ksepg_create_relaxation	(	double	time_delta,
		double	T1,
		double	T2
	)

Create a relaxation operator

Parameters

[in]	time_delta	Time interval over whivh relaxation occurs
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest

Returns

KSEPG_RELAXATION

                                                                                  {

  KSEPG_RELAXATION rel_op;
  rel_op.decayT1 = T1 <= 0 ? 1 : exp(-time_delta/T1);
  rel_op.decayT2 = T2 <= 0 ? 1 : exp(-time_delta/T2);

  return rel_op;

} /* ksepg_create_relaxation() */

ksepg_create_state()

KSEPG_STATE ksepg_create_state

(

int

size

)

Create a new EPG state with a certain maximum number of dephasing steps

Parameters

[in]

size

Maximum number of dephasing steps to be modeled

Returns

KSEPG_STATE

                                         {

  KSEPG_STATE state;
  state.size = size;
  state.data = (double*)calloc(3*size, sizeof(double));

  return state;

} /* ksepg_create_state() */

ksepg_destruct_state()

void ksepg_destruct_state

(

KSEPG_STATE *

state

)

Free an EPG state's memory

Parameters

[in,out]

state

                                              {

  free(state->data);

} /* ksepg_destruct_state() */

ksepg_duplicate_state()

KSEPG_STATE ksepg_duplicate_state

(

const KSEPG_STATE *

state

)

Duplicate an EPG's state

Parameters

[in]

state

to duplicate

Returns

duplicate state

                                                            {

  KSEPG_STATE duplicate = ksepg_create_state(state->size);

  int i = 0;
  for (; i<3*duplicate.size; ++i) {
    duplicate.data[i] = state->data[i];
  }

  return duplicate;

} /* ksepg_duplicate_state() */

ksepg_dephase()

void ksepg_dephase	(	KSEPG_STATE *	state,
		KSEPG_RELAXATION	rel_op
	)

Apply a dephasing step and a relaxation operator to an EPG state

The dephasing step shifts F_+ to the right and F_- to the left.

Parameters

[in,out]	state
[in]	rel_op

                                                                {

  int i;
  int size = state->size;
  double *data = state->data;

  /* shift F_+ to the right */
  for (i=1; i<size; ++i) {
    data[3*(size-i)] = rel_op.decayT2 * data[3*(size-i-1)];
  }

  /* shift F_- to the left */
  for (i=0; i<size-1; ++i) {
    data[3*i + 1] = rel_op.decayT2 * data[3*(i+1) + 1];
  }
  data[3*(size-1) + 1] = 0;

  /* F_+(0) is equal to F_-(0) */
  data[0] = data[1];

  /* apply T1 decay/recovery */
  for (i=0; i<size; ++i) {
    data[3*i + 2] *= rel_op.decayT1;
  }
  data[2] += 1-rel_op.decayT1;

} /* ksepg_dephase() */

ksepg_compute_rotation()

void ksepg_compute_rotation	(	KS_MAT3x3 *	rot,
		double	flip
	)

Compute the rotation matrix corresponding to a refocusing pulse

Parameters

[out]	rot	The computed rotation matrix
[in]	flip	Flip angle of the refocusing pulse

                                                         {

  double c = cos(flip);
  double s = sin(flip);
  double cc = (1+c)/2;
  double ss = (1-c)/2;

  double *r = *rot;
  r[0] = cc;
  r[1] = ss;
  r[2] = -s;

  r[3] = ss;
  r[4] = cc;
  r[5] = s;

  r[6] = s/2;
  r[7] = -s/2;
  r[8] = c;

} /* ksepg_compute_rotation() */

ksepg_apply_rotation()

void ksepg_apply_rotation	(	KSEPG_STATE *	state,
		KS_MAT3x3 *	rot
	)

Apply a rotation matrix to an EPG state

Parameters

[in,out]	state
[in]	rot

                                                                {

  int i, k;
  int size = state->size;
  double *data = state->data;
  const double *rot = *rot_p;
  double tmp[3];

  for (i=0; i<size; ++i) {

    for (k=0; k<3; ++k) {
      tmp[k] = data[3*i]*rot[3*k] + data[3*i+1]*rot[3*k+1] + data[3*i+2]*rot[3*k+2];
    }

    for (k=0; k<3; ++k) {
      data[3*i+k] = tmp[k];
    }
  }

} /* ksepg_apply_rotation() */

ksepg_refocus()

void ksepg_refocus	(	KSEPG_STATE *	state,
		KS_MAT3x3 *	rot
	)

Apply a refocusing operator to an EPG state

This operation differ from ksepg_apply_rotation in that it includes the effect of the crushers.

Parameters

[in,out]	state
[in]	rot	The rotation matrix corresponding to the refocusing pulse

                                                       {

  ksepg_apply_rotation(state, rot);

  /* apply crushers */
  state->data[0] = 0;
  state->data[1] = 0;

} /* ksepg_refocus() */

ksepg_next_echo()

void ksepg_next_echo	(	KSEPG_STATE *	state,
		KS_MAT3x3 *	rot,
		KSEPG_RELAXATION	rel_op
	)

Advance the EPG state from an echo to the next

Parameters

[in,out]	state
[in]	rot	The rotation matrix corresponding to the refocusing pulse between the two echoes
[in]	rel_op	Relaxation operator corresponding to a periof equal to half of the echo spacing

                                                                                  {

  ksepg_dephase(state, rel_op);
  ksepg_refocus(state, rot);
  ksepg_dephase(state, rel_op);

} /* ksepg_next_echo() */

ksepg_get_intensities()

STATUS ksepg_get_intensities	(	const double *	flip_angles,
		int	num_echoes,
		KSEPG_STATE	state,
		double *	signal_intensities,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	state	ADDTEXTHERE
	signal_intensities	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                        {

  if (!flip_angles || !state.data) { return KS_THROW("No flip-angle or state array found"); }

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  int i;
  for (i=0; i<num_echoes; ++i) {
    /* advance the state to before the refocusing */
    ksepg_dephase(&state, rel);

    double flip = flip_angles[i];
    /* ensuring it is not NaN */
    if (isNaN(flip)) {
      return KS_THROW("flip_angles[%d] is NaN", i);
    }

    /* advance the state to the next echo */
    ksepg_compute_rotation(&rot, flip);
    ksepg_refocus(&state, &rot);
    ksepg_dephase(&state, rel);

    /* save refocused signal intensity (state.data[0] and state.data[1] should be equal) */
    if (signal_intensities) {
      signal_intensities[i] = state.data[1];
    }
  }

  return SUCCESS;

} /* ksepg_get_intensities() */

ksepg_get_intensities_std()

STATUS ksepg_get_intensities_std	(	const double *	flip_angles,
		double *	signal_intensities,
		int	num_echoes,
		double	echo_spacing,
		double	T1,
		double	T2
	)

Simulate a classic FSE train with EPG to compute the signal intensities

Parameters

[in]	flip_angles	Refocusing flip angles
[out]	signal_intensities	The intensity levels for each echo relative to T=0
[in]	num_echoes	Number of echoes in the train
[in]	echo_spacing	Echo spacing for the FSE train
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest

Return values

STATUS

SUCCESS or FAILURE

                                                                                            {

  if (!flip_angles || !signal_intensities) { KS_RAISE(FAILURE); }

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* calc signal intensites */
  STATUS status = ksepg_get_intensities(flip_angles, num_echoes,
                                        state,
                                        signal_intensities, 
                                        echo_spacing, T1, T2);

  /* free memory */
  ksepg_destruct_state(&state);

  KS_RAISE(status);

  return SUCCESS;

} /* ksepg_get_intensites_std() */

ksepg_get_intensities_long()

STATUS ksepg_get_intensities_long	(	const double *	flip_angles,
		double *	signal_intensities,
		int	num_echoes,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
	signal_intensities	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                             {

  if (!flip_angles || !signal_intensities) { KS_RAISE(FAILURE); }

  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  /* calc signal intensites */
  STATUS status = ksepg_get_intensities(flip_angles, num_echoes,
                                        state,
                                        signal_intensities, 
                                        echo_spacing, T1, T2);

  /* free memory */
  ksepg_destruct_state(&state);

  KS_RAISE(status);

  return SUCCESS;

} /* ksepg_get_intensities_long() */

ksepg_solve_trig_eq()

double ksepg_solve_trig_eq	(	double	a,
		double	b,
		double	c
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

[in]	a	ADDTEXTHERE
[in]	b	ADDTEXTHERE
	c	ADDTEXTHERE

Return values

value

ADDTEXTHERE

                                                         {

  double R = sqrt(a*a + b*b);
  double beta = atan2(b, a);
  return asin(c/R) - beta;

} /* ksepg_solve_trig_eq() */

ksepg_track_step()

double ksepg_track_step	(	double	Fp,
		double	Fm,
		double	Z,
		double	target_Fm,
		KSEPG_RELAXATION	rel_op
	)

Attempt to track a target intensity

Attempt to select a refocusing flip angle so that F_-(0) will match a desired target intensity at the next echo.

Parameters

[in]	Fp	F_+(1) just before the refocusing pulse
[in]	Fm	F_-(1) just before the refocusing pulse
[in]	Z	Z(1) just before the refocusing pulse
[in]	target_Fm	Target for F_-(0) after the refocusing pulse and a dephasing step
[in]	rel_op	Relaxation operator over the dephasing time

Returns

double Flip angle that tracks the target, NaN if it is not feasible

                                                                                                   {

  double a = Z;
  double b = (Fm - Fp)/2;
  double c = target_Fm/rel_op.decayT2 - (Fm + Fp)/2;

  return ksepg_solve_trig_eq(a, b, c);

} /* ksepg_track_step() */

ksepg_track_intensity_curve()

STATUS ksepg_track_intensity_curve	(	double *	flip_angles,
		int	num_echoes,
		KSEPG_STATE	state,
		double	scaling_factor,
		const double *	target_curve,
		double	echo_spacing,
		double	T1,
		double	T2
	)

Track the requested intensity curve with a custom starting EPG state

This function attempt to set the refocusing flip angle so that the signal intensity of a certain tissue follows over the echoes the requested curve scaled by a factor.

Parameters

[out]	flip_angles	Refocusing flip angle
[in]	num_echoes	Number of echoes in the train
[in,out]	state	EPG state that will be updated
[in]	scaling_factor	Scaling factor for the intensity curve
[in]	target_curve	Target intensity curve to be matched. If NULL it is assumed to be a constant curve equal to one.
[in]	echo_spacing	Echo spacing for the FSE train
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest

Return values

STATUS

SUCCESS or FAILURE

                                                                              {

  if (!flip_angles || !state.data) { return KS_THROW("No flip-angle or state array found"); }

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  int i;
  for (i=0; i<num_echoes; ++i) {
    /* advance the state to before the refocusing */
    ksepg_dephase(&state, rel);

    double flip;
    if (target_curve == NULL) {
      flip = ksepg_track_step(state.data[3], state.data[4], state.data[5], scaling_factor, rel);
    } else {
      flip = ksepg_track_step(state.data[3], state.data[4], state.data[5], target_curve[i]*scaling_factor, rel);
    }

    flip_angles[i] = flip;
    /* ensuring it is not NaN */
    if (isNaN(flip)) {
      return FAILURE;
    }

    /* advance the state to the next echo */
    ksepg_compute_rotation(&rot, flip);
    ksepg_refocus(&state, &rot);
    ksepg_dephase(&state, rel);
  }

  return SUCCESS;

} /* ksepg_track_intensity_curve() */

ksepg_track_intensity_std()

STATUS ksepg_track_intensity_std	(	double *	flip_angles,
		int	num_echoes,
		double	target_intensity,
		double	echo_spacing,
		double	T1,
		double	T2
	)

Track the requested intensity level with a classic FSE train

This function attempt to set the refocusing flip angle so that the signal intensity of a certain tissue is constant over all echoes.

Parameters

[out]	flip_angles	Refocusing flip angle
[in]	num_echoes	Number of echoes in the train
[in]	target_intensity	Target intensity relative to T=0, must be in (0,1)
[in]	echo_spacing	Echo spacing for the FSE train
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest

Return values

STATUS

SUCCESS or FAILURE

                                                                            {

  if (!flip_angles) { return FAILURE; }

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  STATUS s = ksepg_track_intensity_curve(flip_angles, num_echoes,
                                         state,
                                         target_intensity,
                                         NULL,
                                         echo_spacing, T1, T2);
  ksepg_destruct_state(&state);

  return s;

} /* ksepg_track_intensity_std() */

ksepg_track_intensity_long()

STATUS ksepg_track_intensity_long	(	double *	flip_angles,
		int	num_echoes,
		double	target_intensity,
		double	echo_spacing,
		double	T1,
		double	T2
	)

Track the requested intensity level with an extended FSE train

Track the requested intensity level with a modified FSE train where an extra echo with a 90 degree refocusing is inserted at the beginng. It is better suited for long echo trains.

NOTE that the extra echo's data is not meant to be measured.

Parameters

[out]	flip_angles	Refocusing flip angle
[in]	num_echoes	Number of echoes in the train
[in]	target_intensity	Target intensity relative to T=0, must be in (0,1)
[in]	echo_spacing	Echo spacing for the FSE train
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest

Return values

STATUS

SUCCESS or FAILURE

                                                                             {

  if (!flip_angles) { return FAILURE; }

  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  STATUS s = ksepg_track_intensity_curve(flip_angles, num_echoes,
                                         state,
                                         target_intensity,
                                         NULL,
                                         echo_spacing, T1, T2);
  ksepg_destruct_state(&state);

  return s;

} /* ksepg_track_intensity_long() */

ks_print_epg_sim()

void ks_print_epg_sim	(	const float *	flip_angles,
		const double *	si_rel,
		const double *	si_norel,
		const double *	effTEs,
		const int	num_echoes,
		const double	echo_spacing,
		const double	T1,
		const double	T2,
		KS_DESCRIPTION	description
	)

                                                                                                               {

  if (ks_plot_filefmt == KS_PLOT_OFF) {
    return;
  }

  FILE* fp;
  char fname[1000]; /* Full filename of the json file */
  char path[250]; /* Directory of json file */
  char cmd[250];

#ifdef PSD_HW
  /* on MR-scanner */
  sprintf(path,  "/usr/g/mrraw/plot/%s/host/", ks_psdname);
#else
  /* in Simulation */
  sprintf(path, "./plot/host/");
#endif

  sprintf(cmd, "mkdir -p %s > /dev/null", path);
  system(cmd);

  sprintf(fname, "%spsdplot_vrfa_%s_%s.json", path, ks_psdname, description);
  fp = fopen(fname, "w");

  fprintf(fp, "{\n");
  fprintf(fp, "\"psdname\": \"%s\",\n", ks_psdname);
  fprintf(fp, "\"sequence\": \"%s\",\n", description);
  fprintf(fp, "\"T1\": %0.3f,\n", T1);
  fprintf(fp, "\"T2\": %0.3f,\n", T2);
  fprintf(fp, "\"esp\": %0.3f,\n", echo_spacing);
  fprintf(fp, "\"etl\": %d,\n", num_echoes);

  int idx;
  fprintf(fp, "\"flip_angles\": [");
  for (idx = 0; idx < num_echoes - 1; idx++) {
    fprintf(fp, "%f, ", flip_angles[idx]); 
  }
  fprintf(fp, "%f],\n",flip_angles[num_echoes-1]);

  fprintf(fp, "\"si_rel\": [");
  for (idx = 0; idx < num_echoes - 1; idx++) {
    fprintf(fp, "%f, ", si_rel[idx]); 
  }
  fprintf(fp, "%f],\n",si_rel[num_echoes-1]);

  fprintf(fp, "\"si_norel\": [");
  for (idx = 0; idx < num_echoes - 1; idx++) {
    fprintf(fp, "%f, ", si_norel[idx]); 
  }
  fprintf(fp, "%f],\n",si_norel[num_echoes-1]);

  fprintf(fp, "\"effTEs\": [");
  for (idx = 0; idx < num_echoes - 1; idx++) {
    fprintf(fp, "%f, ", effTEs[idx]); 
  }
  fprintf(fp, "%f]\n",effTEs[num_echoes-1]);
  fprintf(fp, "}\n");

  fflush(fp);
  fclose(fp);

} /* ks_print_epg_sim() */

ks_predownload_plot_epg_sim()

void ks_predownload_plot_epg_sim

(

KS_DESCRIPTION

description

)

ADDTITLEHERE

ADDDESCHERE

Parameters

[in]

description

ADDTEXTHERE

Returns

void

                                                             {

  if (ks_plot_filefmt == KS_PLOT_OFF) {
    return;
  }

  char fname[1000]; /* Full filename of the json file */
  char path[250]; /* Directory of json file */
  char outputdir[250]; /* Where plots are saved */
  char outputdir_uid[250]; /* rhkacq_uid based output directory (for automatic transfers) */

#ifdef PSD_HW
  /* on MR-scanner */
  sprintf(path,  "/usr/g/mrraw/plot/%s/host/", ks_psdname);
  sprintf(outputdir, "/usr/g/mrraw/plot/%s/", ks_psdname); 

  if (ks_plot_kstmp) {
    sprintf(outputdir_uid, "/usr/g/mrraw/kstmp/%010d/plot/", rhkacq_uid);
  }
#else
  /* in Simulation */
  sprintf(path, "./plot/host/");
  sprintf(outputdir, "./plot/"); 
#endif

  sprintf(fname, "%spsdplot_vrfa_%s_%s.json", path, ks_psdname, description);

  char cmd[250];

#ifdef PSD_HW
  /* on MR-scanner */
  if (ks_plot_kstmp) {
    sprintf(cmd, "mkdir -p %s > /dev/null", outputdir_uid);
    system(cmd);
  }
  char pythonpath[] = "/usr/local/bin/apython";
  char psdplotpath[] = "/usr/local/bin/psdplot_vrfa.py";
#else
  /* in Simulation */
  char pythonpath[] = "apython";
  char psdplotpath[] = "../KSFoundation/psdplot/psdplot_vrfa.py";
#endif
  sprintf(cmd, "mkdir -p %s > /dev/null", path);
  system(cmd);

  {
  char cmd[1000];
  char fname_plan[1000]; /* Full filename of the json file */
  sprintf(fname_plan, "%spsdplot_%s_%s.json", path, ks_psdname, description);
  sprintf(cmd, "%s %s %s %s "
                "--outputdir %s &",
          pythonpath, psdplotpath, fname_plan, fname,
          outputdir);
  ks_dbg("%s", cmd);
  system(cmd);
#ifdef PSD_HW
  /* on MR-scanner */
  if (ks_plot_kstmp) {
    sprintf(outputdir_uid, "/usr/g/mrraw/kstmp/%010d/plot/", rhkacq_uid);
    sprintf(cmd, "cp %s%s_%s.html %s",
          outputdir, ks_psdname, description,
          outputdir_uid);
    system(cmd);
    sprintf(cmd, "cp %s %s > /dev/null", fname, outputdir_uid);
    system(cmd);
  }
#endif
  }
(void)outputdir_uid;

} /* ks_predownload_plot_epg_sim() */

ksepg_calc_effTE()

STATUS ksepg_calc_effTE	(	const float *	flip_angles,
		double *	effTEs,
		int	num_echoes,
		double	echo_spacing,
		double	T1,
		double	T2,
		KS_DESCRIPTION	desc
	)

Computes the effective echo time at each echo of a classic FSE train

Parameters

[in]	flip_angles	Refocusing flip angles
[out]	effTEs	Effective echo time for each echo
[in]	num_echoes	Number of echoes in the train
[in]	echo_spacing	Echo spacing for the FSE train
[in]	T1	T_1 value for the tissue of interest
[in]	T2	T_2 value for the tissue of interest
[in]	desc	Description for printing to file

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                                         {

  if (!flip_angles || !effTEs) { KS_RAISE(FAILURE); }

  STATUS status;
  int i; 

  double *flip_angles_radians = (double *)alloca(num_echoes * sizeof(double));
  for (i=0; i<num_echoes; i++){
    flip_angles_radians[i] = (double)flip_angles[i]*PI/180;
  }

  /* calc signal intensites with and without relaxation */
  double *si_rel = (double *)alloca(num_echoes * sizeof(double));
  status = ksepg_get_intensities_std(flip_angles_radians, si_rel, num_echoes, echo_spacing, T1, T2);
  KS_RAISE(status);

  ks_dbg("intensity (esp=%.1f, T1=%.0f, T2=%.0f): %.2f -> %.2f -> %.2f",
         echo_spacing, T1, T2, si_rel[1], si_rel[num_echoes/2], si_rel[num_echoes-1]);

  double *si_norel = (double *)alloca(num_echoes * sizeof(double));
  status = ksepg_get_intensities_std(flip_angles_radians, si_norel, num_echoes, 0, 1, 1);
  KS_RAISE(status);

  /* calculate relaxation term (r) and effTEs */
  int echo; 
  for (echo=0; echo < num_echoes; echo++) {

    if (si_rel[echo]>0 && si_norel[echo]>=si_rel[echo]) {

      double r = si_rel[echo] / si_norel[echo];
      effTEs[echo] = -T2*log(r);
    } else {

      /* Failed to calculate, use a placeholder instead */
      effTEs[echo] = (echo+1)*echo_spacing;
    }
  }

  ks_print_epg_sim(flip_angles, si_rel, si_norel, effTEs, num_echoes, echo_spacing, T1, T2, description);

  /* Check for "slice profile" */
  const double rel_beg = si_rel[1];
  const double rel_mid = si_rel[num_echoes/2];
  const double rel_end = si_rel[num_echoes-1];


  for (i=0; i<num_echoes; i++){
    flip_angles_radians[i] *= 0.5;
  }

  status = ksepg_get_intensities_std(flip_angles_radians, si_rel, num_echoes, 0, 1, 1);
  KS_RAISE(status);
  ks_dbg("intensity half norel loss: %.2f -> %.2f -> %.2f",
         si_rel[1]/si_norel[1], si_rel[num_echoes/2]/si_norel[num_echoes/2], si_rel[num_echoes-1]/si_norel[num_echoes-1]);

  status = ksepg_get_intensities_std(flip_angles_radians, si_rel, num_echoes, echo_spacing, T1, T2);
  KS_RAISE(status);

  ks_dbg("intensity half rel loss:   %.2f -> %.2f -> %.2f",
         si_rel[1]/rel_beg, si_rel[num_echoes/2]/rel_mid, si_rel[num_echoes-1]/rel_end);  

  return SUCCESS;

} /* ksepg_calc_effTE() */

ksepg_max_angle()

double ksepg_max_angle	(	double *	flip_angles,
		int	num_echoes
	)

Returns the maximum refocusing angle

Parameters

[in]	flip_angles	Refocusing flip angles
[in]	num_echoes	Number of echoes in the train

Returns

double The maximum refocusing flip angle

                                                            {

  if (!flip_angles) { return 0; }

  double max_angle = 0;
  int i;
  for (i=0; i<num_echoes; ++i) {
    max_angle = flip_angles[i] > max_angle ? flip_angles[i] : max_angle;
  }

  return max_angle;

} /* ksepg_max_angle() */

ksepg_optimize_constant()

STATUS ksepg_optimize_constant	(	double *	intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		const KSEPG_STATE	starting_state
	)

ADDTITLEHERE

Parameters

	intensity	ADDTEXTHERE
	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE
[in]	starting_state	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                 {

  if (!flip_angles) { return KS_THROW("NULL input"); }

  double low_tgt = 0;
  double hi_tgt = 1;

  STATUS status = FAILURE;
  int i;
  for (i=0; i<10; ++i) {
    double tgt = (hi_tgt + low_tgt)/2;

    KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
    status = ksepg_track_intensity_curve(flip_angles, num_echoes, state, tgt, NULL, echo_spacing, T1, T2);
    ksepg_destruct_state(&state);

    if (status == SUCCESS && ksepg_max_angle(flip_angles, num_echoes) < max_flip) {
      low_tgt = tgt;
    } else {
      hi_tgt = tgt;
    }
  }

  if (low_tgt > 0) {
    if (intensity) {
      *intensity = low_tgt;
    }

    KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
    status = ksepg_track_intensity_curve(flip_angles, num_echoes, state, low_tgt, NULL, echo_spacing, T1, T2);
    ksepg_destruct_state(&state);
    KS_RAISE(status);

    return SUCCESS;
  }

  return KS_THROW("Failed");

} /* ksepg_optimize_constant() */

ksepg_optimize_constant_std()

STATUS ksepg_optimize_constant_std	(	double *	intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2
	)

Find the maximum constant intensity that can be mantained

                                                                              {

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  return ksepg_optimize_constant(intensity,
                                 flip_angles, num_echoes,
                                 max_flip,
                                 echo_spacing, T1, T2,
                                 state);

} /* ksepg_optimize_constant_std() */

ksepg_optimize_constant_long()

STATUS ksepg_optimize_constant_long	(	double *	intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	intensity	ADDTEXTHERE
	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                               {

  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  return ksepg_optimize_constant(intensity,
                                 flip_angles, num_echoes,
                                 max_flip,
                                 echo_spacing, T1, T2,
                                 state);

} /* ksepg_optimize_constant_long() */

ksepg_generate_intensity_triangle_impl()

STATUS ksepg_generate_intensity_triangle_impl	(	double *	target_intensities,
		int	num_echoes,
		double	center_position,
		double	min_intensity,
		double	max_intensity,
		double	half_size
	)

Circular VFA modulation

The aim is to modulate the FA so that the intensity is the same at the start and the end of the train. There are two possible strategies:

• Slanted triangle, where the peak is where the center echo is intented to be and both sides decrease with different slopes towards the same intensity.
• Circular shifted triangle, where a symmetric triangle is circulally shifted in order to match the peak with the center echo. There should be limitations on where the center echo can be dispaced from the middle. Maybe displacements greater than num_echoes/4 should not be allowed. The intensity at the center echoes is an input. The intensity at the valleys of the triangle is optimized. In general we want the highest intensity possible up to matching the center echo's

Parameters

	target_intensities	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	center_position	ADDTEXTHERE
[in]	min_intensity	ADDTEXTHERE
[in]	max_intensity	ADDTEXTHERE
[in]	half_size	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                                    {
  if (!target_intensities) {
    return KS_THROW("Target_intensities is NULL");
  }

  if (center_position < 0 || center_position >= num_echoes) {
    return KS_THROW("center_echo (%f) must be in the interval [0, %d]", center_position, num_echoes-1);
  }

  if (min_intensity < 0 || min_intensity >= 1) {
    return KS_THROW("min_intensity (%f) must be in the interval [0, 1)", min_intensity);
  }

  if (max_intensity <= min_intensity || max_intensity >= 1) {
    return KS_THROW("max_intensity (%f) must be in the interval [%f, 1)", max_intensity, min_intensity);
  }

  const double min_half_size = (num_echoes - 1.0)/2.0;
  if (half_size < min_half_size) {
    return KS_THROW("half_size (%f) must be greater than half of num_echoes %f", half_size, min_half_size);
  }

  int i;

  /* Generate the pattern */
  double max_unscaled = 0, min_unscaled = 1;
  for (i=0; i<num_echoes; ++i) {
    double echo_pos = i - center_position;
    double unscaled = fabs(1 - fabs(echo_pos)/half_size);

    max_unscaled = FMax(2, unscaled, max_unscaled);
    min_unscaled = FMin(2, unscaled, min_unscaled);
    target_intensities[i] = unscaled;
  }

  /* Scale the pattern to min_intensity/max_intensity */
  for (i=0; i<num_echoes; ++i) {
    target_intensities[i] = min_intensity + (target_intensities[i] - min_unscaled)*(max_intensity - min_intensity) / (max_unscaled - min_unscaled);
  }

  return SUCCESS;

} /* ksepg_generate_intensity_triangle_impl() */

ksepg_generate_intensity_periodic_triangle()

STATUS ksepg_generate_intensity_periodic_triangle	(	double *	target_intensities,
		int	num_echoes,
		double	center_position,
		double	min_intensity,
		double	max_intensity
	)

Function used to generate bounded intensity curves

                                                                                                                      {
  const double half_size = (num_echoes - 1.0)/2.0;

  return ksepg_generate_intensity_triangle_impl(target_intensities, num_echoes,
                                                center_position, min_intensity, max_intensity, half_size);

} /* ksepg_generate_intensity_periodic_triangle() */

ksepg_generate_intensity_cropped_triangle()

STATUS ksepg_generate_intensity_cropped_triangle	(	double *	target_intensities,
		int	num_echoes,
		double	center_position,
		double	min_intensity,
		double	max_intensity
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	target_intensities	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	center_position	ADDTEXTHERE
[in]	min_intensity	ADDTEXTHERE
[in]	max_intensity	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                                     {
  const double half_size = FMax(2, center_position, num_echoes - 1 - center_position);

  return ksepg_generate_intensity_triangle_impl(target_intensities, num_echoes,
                                                center_position, min_intensity, max_intensity, half_size);

} /* ksepg_generate_intensity_cropped_triangle() */

ksepg_generate_intensity_ramp()

STATUS ksepg_generate_intensity_ramp	(	double *	target_intensities,
		int	num_echoes,
		double	center_position,
		double	min_intensity,
		double	max_intensity
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	target_intensities	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	center_position	ADDTEXTHERE
[in]	min_intensity	ADDTEXTHERE
[in]	max_intensity	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                         {

  int i = 0;
  for (; i < center_position+1; ++i) {
    target_intensities[i] = max_intensity;
  }
  const double hedge = i-1;
  for (; i<num_echoes; ++i) {
    target_intensities[i] =
      max_intensity * (num_echoes-1 - i)/(num_echoes-1 - hedge) +
      min_intensity * (i - hedge)/(num_echoes-1 - hedge);
  }

  return SUCCESS;

} /* ksepg_generate_intensity_ramp() */

ksepg_generate_intensity_linear()

STATUS ksepg_generate_intensity_linear	(	double *	target_intensities,
		int	num_echoes,
		double	center_position,
		double	min_intensity,
		double	max_intensity
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	target_intensities	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	center_position	ADDTEXTHERE
[in]	min_intensity	ADDTEXTHERE
[in]	max_intensity	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                                           {

  (void)center_position;
  double slope = (min_intensity-max_intensity)/(num_echoes-1);
  int i = 0;
  for (; i<num_echoes; ++i) {
    target_intensities[i] = max_intensity + i*slope;
  }

  return SUCCESS;

} /* ksepg_generate_intensity_linear() */

ksepg_track_bounded_curve()

STATUS ksepg_track_bounded_curve	(	double *	flip_angles,
		double *	target_intensities,
		int	num_echoes,
		double	min_intensity,
		double	max_intensity,
		double	echo_spacing,
		double	T1,
		double	T2,
		double	center_position,
		STATUS(*)(double *, int, double, double, double)	generate_bounded_curve,
		const KSEPG_STATE	starting_state
	)

                                                                   {
  STATUS s;
  s = generate_bounded_curve(target_intensities, num_echoes, center_position, min_intensity, max_intensity);
  KS_RAISE(s);

  KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
  s = ksepg_track_intensity_curve(flip_angles, num_echoes, state,
                                  1.0, target_intensities,
                                  echo_spacing, T1, T2);
  ksepg_destruct_state(&state);

  return s;

} /* ksepg_track_bounded_curve() */

ksepg_optimize_bounded_curve()

STATUS ksepg_optimize_bounded_curve	(	double *	min_intensity,
		double *	max_intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		double	center_position,
		STATUS(*)(double *, int, double, double, double)	generate_bounded_curve,
		const KSEPG_STATE	starting_state
	)

                                                                      {

  if (!flip_angles) { return KS_THROW("No flip-angle or state array found"); }

  STATUS status = FAILURE;
  double opt_max_intensity = 0;
  double opt_min_intensity = 0;

  double *target_intensities;
  target_intensities = (double*)alloca(num_echoes * sizeof(double));
  double opt_objective = 0.0;

  double alpha; 
  const double outer_loop_step = 0.1;
  int i;
  const int inner_loop_size = 20;

  if (outer_loop_step < 0.01 || outer_loop_step > 0.5) {
    return KS_THROW("outer_loop_step (%f) is out of bounds", outer_loop_step);
  }

  for (alpha=outer_loop_step; alpha<1.0; alpha += outer_loop_step) {

    double max_high = 1;
    double max_low = 0;

    for (i=0; i<inner_loop_size; i++) {

      const double loop_max_intensity = (max_high + max_low)/2;
      const double loop_min_intensity = alpha * loop_max_intensity;

      status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                         loop_min_intensity, loop_max_intensity,
                                         echo_spacing, T1, T2,
                                         center_position,
                                         generate_bounded_curve,
                                         starting_state);

      if (status != SUCCESS  || ksepg_max_angle(flip_angles, num_echoes) > max_flip) {
        max_high = loop_max_intensity;
        continue;
      }

      max_low = loop_max_intensity;

      const double objective = pow(loop_max_intensity, 4) * loop_min_intensity;
      if (opt_objective >= objective) {
        continue;
      }

      opt_objective = objective;
      opt_max_intensity = loop_max_intensity;
      opt_min_intensity = loop_min_intensity;
    }
  }

  status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                     opt_min_intensity, opt_max_intensity,
                                     echo_spacing, T1, T2,
                                     center_position,
                                     generate_bounded_curve,
                                     starting_state);
  KS_RAISE(status);

  if (min_intensity) { *min_intensity = opt_min_intensity; };
  if (max_intensity) { *max_intensity = opt_max_intensity; };

  return SUCCESS;

} /* ksepg_optimize_bounded_curve */

ksepg_optimize_bounded_curve_std()

STATUS ksepg_optimize_bounded_curve_std	(	double *	min_intensity,
		double *	max_intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		double	center_position,
		STATUS(*)(double *, int, double, double, double)	generate_bounded_curve
	)

Optimize a bounded intensity curve so that it maximizes the center position's intensity

Optimize the target intensity curve modeled as a bounded curve in order to maximise the intensity at the echo corresponding to the center of the k-space. A bounded curve is parameterized by the maximum and minim intensities over all echoes.

                                                                                                                {
  const int state_size = 2*num_echoes + 1;
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  STATUS s = ksepg_optimize_bounded_curve(min_intensity, max_intensity,
                                          flip_angles, num_echoes,
                                          max_flip,
                                          echo_spacing, T1, T2,
                                          center_position,
                                          generate_bounded_curve,
                                          state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_bounded_curve_std() */

ksepg_optimize_bounded_curve_long()

STATUS ksepg_optimize_bounded_curve_long	(	double *	min_intensity,
		double *	max_intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		double	center_position,
		STATUS(*)(double *, int, double, double, double)	generate_bounded_curve
	)

ADDTITLEHERE

                                                                                                                 {
  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  STATUS s = ksepg_optimize_bounded_curve(min_intensity, max_intensity,
                                          flip_angles, num_echoes,
                                          max_flip,
                                          echo_spacing, T1, T2,
                                          center_position,
                                          generate_bounded_curve,
                                          state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_bounded_curve_long() */

ksepg_optimize_rampdown()

STATUS ksepg_optimize_rampdown	(	double *	end_intensity,
		double	start_intensity,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		const KSEPG_STATE	starting_state
	)

ADDTITLEHERE

                                                                 {

  if (!flip_angles) { return KS_THROW("No flip-angle or state array found"); }

  STATUS status = FAILURE;
  *end_intensity = 1;

  double *target_intensities;
  target_intensities = (double *)alloca(num_echoes * sizeof(double));

  int i;
  const int loop_size = 10;

  double end_high = 1;
  double end_low = 0;
  for (i=0; i<loop_size; i++) {

    const double loop_end_intensity = (end_high + end_low)/2;

    status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                       loop_end_intensity, start_intensity,
                                       echo_spacing, T1, T2,
                                       0,
                                       ksepg_generate_intensity_linear,
                                       starting_state);

    if (status != SUCCESS  || ksepg_max_angle(flip_angles, num_echoes) > max_flip) {
      end_high = loop_end_intensity;
      continue;
    }

    end_low = loop_end_intensity;

    *end_intensity = loop_end_intensity;
  }

  status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                     *end_intensity, start_intensity,
                                     echo_spacing, T1, T2,
                                     0,
                                     ksepg_generate_intensity_linear,
                                     starting_state);
  KS_RAISE(status);

  return SUCCESS;

} /* ksepg_optimize_rampdown() */

ksepg_half_intensity()

STATUS ksepg_half_intensity	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		KSEPG_TAIL	tail_type
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE
[in]	tail_type	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                  {

  if (!flip_angles) { return KS_THROW("NULL flip-angle array"); }

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  const double intensity = 0.5;

  int i;
  for (i=0; i<num_echoes; ++i) {
    /* advance the state to before the refocusing */
    ksepg_dephase(&state, rel);

    const double flip = ksepg_track_step(state.data[3], state.data[4], state.data[5], intensity, rel);

    if (isNaN(flip) || flip > max_flip) {
      break;
    }

    flip_angles[i] = flip;

    /* advance the state to the next echo */
    ksepg_compute_rotation(&rot, flip);
    ksepg_refocus(&state, &rot);
    ksepg_dephase(&state, rel);
  }

  switch (tail_type) {
  case KSEPG_RAMPDOWN:
    {
      return KS_THROW("Not implemented yet");
    }
  default:
    {
      for (; i<num_echoes; ++i) {
        flip_angles[i] = max_flip;
      }
    }
  }

  return SUCCESS;

} /* ksepg_half_intensity() */

ksepg_optimize_expdown()

STATUS ksepg_optimize_expdown	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	echo_spacing,
		double	T1,
		double	T2,
		const KSEPG_STATE	starting_state
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE
[in]	starting_state	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                {

  if (!flip_angles) { return KS_THROW("NULL flip-angle array"); }

  STATUS status;
  const double start_intensity = starting_state.data[1];

  double alpha_high = 0.1;
  double alpha_low = -0.1;

  double *target_intensity = (double*)alloca(num_echoes * sizeof(double));

  double alpha = alpha_low;
  int i = 0;
  int echo;
  for (; i<10; ++i) {

    /* Generate the target intensity curve */
    for (echo = 0; echo<num_echoes; ++echo) {
      target_intensity[echo] = start_intensity*exp(alpha*(echo + 0.5)*echo_spacing);
    }

    KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
    status = ksepg_track_intensity_curve(flip_angles, num_echoes,
                                         state,
                                         1, target_intensity,
                                         echo_spacing, T1, T2);
    if (status != SUCCESS && ksepg_max_angle(flip_angles, num_echoes) < max_flip) {
      alpha_high = alpha;
      continue;
    }

    alpha_low = alpha;
    alpha = (alpha_high + alpha_low)/2;
  }

  /* Rerun the EPG to fill the correct flip angles */
  for (echo=0; echo<num_echoes; ++echo) {
    target_intensity[echo] = start_intensity*exp(alpha_low*(echo + 0.5)*echo_spacing);
  }

  KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
  status = ksepg_track_intensity_curve(flip_angles, num_echoes,
                                       state,
                                       1, target_intensity,
                                       echo_spacing, T1, T2);
  KS_RAISE(status);

  return SUCCESS;

} /* ksepg_optimize_expdown() */

ksepg_optimize_half_exponential()

STATUS ksepg_optimize_half_exponential	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                        {

  if (!flip_angles) { return KS_THROW("NULL flip-angle array"); }

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(0, 1, 1);

  /* apply a first refocusing of 90 degrees */
  flip_angles[0] = PI/2;
  ksepg_compute_rotation(&rot, flip_angles[0]);
  ksepg_next_echo(&state, &rot, rel);

  STATUS s = ksepg_optimize_expdown(flip_angles+1, num_echoes-1,
                                    max_flip,
                                    0, 1, 1,
                                    state);

  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_half_exponential() */

ksepg_optimize_half_free()

STATUS ksepg_optimize_half_free	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                 {

  if (!flip_angles) { return KS_THROW("NULL flip-angle array"); }

  const int state_size = 2*num_echoes + 1;

  /* initialize the state to the value after the excitation */
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(0, 1, 1);

  /* apply a first refocusing of 90 degrees */
  flip_angles[0] = PI/2;
  ksepg_compute_rotation(&rot, flip_angles[0]);
  ksepg_next_echo(&state, &rot, rel);

  const double start_intensity = state.data[1];
  double end_intensity;

  STATUS s = ksepg_optimize_rampdown(&end_intensity, start_intensity,
                                     flip_angles+1, num_echoes-1,
                                     max_flip,
                                     0, 1, 1,
                                     state);

  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_half_free() */

ksepg_optimize_linear_free()

STATUS ksepg_optimize_linear_free	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2,
		const KSEPG_STATE	starting_state
	)

Search for a linear relaxation-free intensity curve that

maximize the concave objective function:

f_obj = max_intensity * min_intensity^min_intensity_weight

where both max_intensity and min_intensity are calculated by taking into accout the relexation effect.

Parameters

[out]	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE
[in]	starting_state	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                    {


  if (!flip_angles) {
    return KS_THROW("No flip-angle or state array found");
  }
  if (min_intensity_weight < 0) {
    return KS_THROW("Weight for the minimum intensity must be non-negative");
  }

  STATUS status = FAILURE;
  double opt_max_intensity = 0;
  double opt_min_intensity = 0;

  double *target_intensities;
  target_intensities = (double*)alloca(num_echoes * sizeof(double));
  double opt_objective = 0.0;

  double max_intensity; 
  const double outer_loop_step = 0.1;
  int i;
  const int inner_loop_size = 10;

  if (outer_loop_step < 0.01 || outer_loop_step > 0.5) {
    return KS_THROW("outer_loop_step (%f) is out of bounds", outer_loop_step);
  }

  for (max_intensity=outer_loop_step; max_intensity<1.0; max_intensity += outer_loop_step) {

    double min_high = 1;
    double min_low = max_intensity;

    for (i=0; i<inner_loop_size; i++) {
      const double min_intensity = (min_high + min_low)/2;

      status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                         min_intensity, max_intensity,
                                         0, 1, 1, /* relaxation free */
                                         0, /* unused */
                                         ksepg_generate_intensity_linear,
                                         starting_state);

      if (status != SUCCESS  || ksepg_max_angle(flip_angles, num_echoes) > max_flip) {
        min_high = min_intensity;
        continue;
      }

      min_low = min_intensity;

      /* Calculate the objective by */
      KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
      status = ksepg_get_intensities(flip_angles, num_echoes,
                                     state,
                                     target_intensities,
                                     echo_spacing, T1, T2);
      ksepg_destruct_state(&state);
      KS_RAISE(status);

      const double objective =
        target_intensities[0] * pow(target_intensities[num_echoes-1], min_intensity_weight);  
      if (opt_objective > objective) {
        continue;
      }

      opt_objective = objective;
      opt_max_intensity = max_intensity;
      opt_min_intensity = min_intensity;
    }
  }

  status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                     opt_min_intensity, opt_max_intensity,
                                     0, 1, 1, /* relaxation free */
                                     0, /* unused */
                                     ksepg_generate_intensity_linear,
                                     starting_state);
  KS_RAISE(status);

  return SUCCESS;

} /* ksepg_optimize_linear_free() */

ksepg_optimize_linear_free_std()

STATUS ksepg_optimize_linear_free_std	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                 {
  const int state_size = 2*num_echoes + 1;
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  STATUS s = ksepg_optimize_linear_free(flip_angles, num_echoes,
                                        max_flip, min_intensity_weight,
                                        echo_spacing, T1, T2,
                                        state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;
  
} /* ksepg_optimize_linear_free_std() */

ksepg_optimize_linear_free_long()

STATUS ksepg_optimize_linear_free_long	(	double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                  {
  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  STATUS s = ksepg_optimize_linear_free(flip_angles, num_echoes,
                                        max_flip, min_intensity_weight,
                                        echo_spacing, T1, T2,
                                        state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_linear_free_long() */

ksepg_optimize_linear_free_center_echo()

STATUS ksepg_optimize_linear_free_center_echo	(	int *	_center_echo,
		int	target_eff_TE,
		int	linear_sweep_center_echo,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2,
		const KSEPG_STATE	starting_state
	)

Joint VFA modulation and center echo selection from TE

The objective is to place the center and modulate the FA so that we get a certain effective TE.

Start placing the center in the middle, target an VFA that would yield an intensity that matches the SE intensity at the requested TE for certain T1/T2 values. Calculate the effective TE. Move the center to the left or right depending if the effective TE is higher or lower than the target TE respectively. Regenerate the VFA modulation with the same target for the center echo.

Does this converge?

Parameters

	_center_echo	ADDTEXTHERE
[in]	target_eff_TE	ADDTEXTHERE
[in]	linear_sweep_center_echo	ADDTEXTHERE
	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE
[in]	starting_state	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                {


  if (!flip_angles) {
    return KS_THROW("No flip-angle or state array found");
  }
  if (min_intensity_weight < 0) {
    return KS_THROW("Weight for the minimum intensity must be non-negative");
  }

  STATUS status = FAILURE;
  double opt_max_intensity = 0;
  double opt_min_intensity = 0;
  int opt_center_echo = 0;

  double *target_intensities;
  target_intensities = (double*)alloca(num_echoes * sizeof(double));
  double *relax_intensities;
  relax_intensities = (double*)alloca(num_echoes * sizeof(double));
  double opt_objective = 0.0;

  double max_intensity; 
  const double outer_loop_step = 0.1;
  int i;
  const int inner_loop_size = 20;

  if (outer_loop_step < 0.01 || outer_loop_step > 0.5) {
    return KS_THROW("outer_loop_step (%f) is out of bounds", outer_loop_step);
  }

  const double normalizing_factor = IMin(2, linear_sweep_center_echo, num_echoes - linear_sweep_center_echo -1);

  for (max_intensity=outer_loop_step; max_intensity<1.0; max_intensity += outer_loop_step) {

    double min_high = 1;
    double min_low = 0;

    for (i=0; i<inner_loop_size; i++) {
      const double min_intensity = (min_high + min_low)/2;

      status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                         min_intensity, max_intensity,
                                         0, 1, 1, /* relaxation free */
                                         0, /* unused */
                                         ksepg_generate_intensity_linear,
                                         starting_state);

      if (status != SUCCESS  || ksepg_max_angle(flip_angles, num_echoes) > max_flip) {
        min_high = min_intensity;
        continue;
      }

      min_low = min_intensity;

      /* Calculate the objective by */
      KSEPG_STATE state = ksepg_duplicate_state(&starting_state);
      status = ksepg_get_intensities(flip_angles, num_echoes,
                                     state,
                                     relax_intensities,
                                     echo_spacing, T1, T2);
      ksepg_destruct_state(&state);
      KS_RAISE(status);

      int echo, center_echo;
      for (echo=0; echo < num_echoes; ++echo) {
        double r = relax_intensities[echo] / target_intensities[echo];
        double eff_TE = -T2 * log(r);
        if (eff_TE > target_eff_TE) {
          center_echo = echo > 1 ? echo -1 : 0;
          break;
        }
      }

      const double relative_distance = (center_echo - linear_sweep_center_echo)/normalizing_factor;
      const double objective = (relax_intensities[0] + 1) *
        (pow(relax_intensities[num_echoes-1], min_intensity_weight) + 1) /
        (pow(relative_distance, 2) + 1);  
      if (opt_objective >= objective) {
        continue;
      }

      opt_objective = objective;
      opt_center_echo = center_echo;
      opt_max_intensity = max_intensity;
      opt_min_intensity = min_intensity;
    }
  }

  status = ksepg_track_bounded_curve(flip_angles, target_intensities, num_echoes,
                                     opt_min_intensity, opt_max_intensity,
                                     0, 1, 1, /* relaxation free */
                                     0, /* unused */
                                     ksepg_generate_intensity_linear,
                                     starting_state);
  KS_RAISE(status);

  if (_center_echo) { *_center_echo = opt_center_echo; };

  return SUCCESS;

} /* ksepg_optimize_linear_free_center_echo() */

ksepg_optimize_linear_free_center_echo_std()

STATUS ksepg_optimize_linear_free_center_echo_std	(	int *	center_echo,
		int	target_eff_TE,
		int	linear_sweep_center_echo,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	center_echo	ADDTEXTHERE
[in]	target_eff_TE	ADDTEXTHERE
[in]	linear_sweep_center_echo	ADDTEXTHERE
[in]	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                             {
  const int state_size = 2*num_echoes + 1;
  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  STATUS s = ksepg_optimize_linear_free_center_echo(center_echo,
                                                    target_eff_TE, linear_sweep_center_echo,
                                                    flip_angles, num_echoes,
                                                    max_flip, min_intensity_weight,
                                                    echo_spacing, T1, T2,
                                                    state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_linear_free_center_echo_std() */

ksepg_optimize_linear_free_center_echo_long()

STATUS ksepg_optimize_linear_free_center_echo_long	(	int *	center_echo,
		int	target_eff_TE,
		int	linear_sweep_center_echo,
		double *	flip_angles,
		int	num_echoes,
		double	max_flip,
		double	min_intensity_weight,
		double	echo_spacing,
		double	T1,
		double	T2
	)

ADDTITLEHERE

ADDDESCHERE

Parameters

	center_echo	ADDTEXTHERE
[in]	target_eff_TE	ADDTEXTHERE
[in]	linear_sweep_center_echo	ADDTEXTHERE
	flip_angles	ADDTEXTHERE
[in]	num_echoes	ADDTEXTHERE
[in]	max_flip	ADDTEXTHERE
[in]	min_intensity_weight	ADDTEXTHERE
[in]	echo_spacing	ADDTEXTHERE
[in]	T1	ADDTEXTHERE
[in]	T2	ADDTEXTHERE

Return values

STATUS

SUCCESS or FAILURE

                                                                                              {
  const int state_size = 2*num_echoes + 3;

  KS_MAT3x3 rot;
  KSEPG_RELAXATION rel = ksepg_create_relaxation(echo_spacing/2, T1, T2);

  KSEPG_STATE state = ksepg_create_state(state_size);
  state.data[0] = 1;
  state.data[1] = 1;

  /* apply a first refocusing of 90 degrees */
  ksepg_compute_rotation(&rot, PI/2);
  ksepg_next_echo(&state, &rot, rel);

  STATUS s = ksepg_optimize_linear_free_center_echo(center_echo,
                                                    target_eff_TE, linear_sweep_center_echo,
                                                    flip_angles, num_echoes,
                                                    max_flip, min_intensity_weight,
                                                    echo_spacing, T1, T2,
                                                    state);
  ksepg_destruct_state(&state);

  KS_RAISE(s);

  return SUCCESS;

} /* ksepg_optimize_linear_free_center_echo_long() */

ksepg_sineTRAPS()

STATUS ksepg_sineTRAPS	(	double *	flip_angles,
		int	etl,
		int	n23,
		int	n4,
		double	an1,
		double	an23,
		double	an4
	)

Generates an array of refocusing flip angles in the shape of a sinusoidal TRAPS scheme

Legacy function, please use ksepg_sineTRAPS2 instead.

From the following paper: Weigel, M., & Hennig, J. (2006). Contrast behavior and relaxation effects of conventional and hyperecho-turbo spin echo sequences at 1.5 and 3 T. Magnetic Resonance in Medicine: Official Journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, 55(4), 826–835.

Parameters

[out]	flip_angles	Array of flip angles [degrees]
[in]	etl	Number of refocusing pulses
[in]	n23	Index at the FA peak
[in]	n4	Index at the start of the tail
[in]	an1	RF amplitude at the start of the first ramp [degrees]
[in]	an23	RF amplitude on the plateau [degrees]
[in]	an4	RF amplitude on the tail [degrees]

Return values

STATUS

SUCCESS or FAILURE

                                   {


  double scale_fact = 0.92;
  int n1 = 2;
  int index = 0;
  int i; 

  /* Calculate length of stuff */
  int tail = etl - 2;
  int n_sine_a = n23 - n1;
  int n_flat = etl - n4;
  int n_sine_d = tail - n_sine_a + 1 - n_flat;

  /* TODO: error checking */

  /* Ramp */
  double a1 = 0.6;
  double a2 = 0.1;

  flip_angles[index++] = an1+a1*(180.0-an1);
  flip_angles[index++] = an1+a2*(180.0-an1);

  /* Acending sine */
  an23 /= scale_fact;
  for (i = 0; i < n_sine_a; i++) {

    double radi = (-PI/2.0) + ((PI/2.0)/(n_sine_a-1))*i;

    double flip_ratio = 1.0 + (((an23/an1/2.0)-1)/(n_sine_a-1))*i;
    flip_ratio *= an1;

    flip_angles[index++] = (sin(radi)+2) * flip_ratio;

  }
  flip_angles[index-1] *= scale_fact;

  /* Decending sine */
  for (i = 0; i < n_sine_d; i++) {

    double radi = -PI + ((PI/2.0)/(n_sine_d-1))*i;

    double start = (an23/an4/2);
    double flip_ratio = start + ((1-start)/(n_sine_d-1))*i;
    flip_ratio *= an4;

    if (i != 0) {
      flip_angles[index++] = (sin(radi)+2) * flip_ratio;
    }

  }

  /* Flat */
  for (i = 0; i < n_flat; i++) {
    flip_angles[index++] = an4;
  }

  return SUCCESS; 

} /* ksepg_sineTRAPS() */

ksepg_sineTRAPS2()

STATUS ksepg_sineTRAPS2	(	float *	flip_angles,
		int	etl,
		KSEPG_TRAPS_DESIGN	traps_design,
		float	scaling_factor
	)

Generates an array of refocusing flip angles in the shape of a sinusoidal TRAPS scheme

From the following paper: Weigel, M., & Hennig, J. (2006). Contrast behavior and relaxation effects of conventional and hyperecho-turbo spin echo sequences at 1.5 and 3 T. Magnetic Resonance in Medicine: Official Journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine, 55(4), 826–835.

Parameters

[out]	flip_angles	Array of flip angles [degrees]
[in]	etl	Number of refocusing pulses
[in]	traps_design	Design structure
[in]	scaling_factor	Scaling factor for all flip angles (1.0 recommended)

Return values

STATUS

SUCCESS or FAILURE

                                             {

  /* rule of thumb the lower an1 the longer the ramp needed to get into pseudo steady state (PSS) */                       
  float ramp[3];
  switch(traps_design.ramp_length) {
    case 1: 
      ramp[0] = 0.45;
      break;
    case 2:
      ramp[0] = 0.6;
      ramp[1] = 0.1;
      break;
    case 3:
      ramp[0] = 0.6;
      ramp[1] = 0.1;
      ramp[2] = 0.05;
      break;
    default: 
      return KS_THROW("ramp_length must be 1,2 or 3");  
  }

  if (traps_design.n1 < traps_design.ramp_length || traps_design.n1 >= traps_design.n2 /*space for n23 and n4*/){
    return KS_THROW("n1 must fullfill ramp_length <= n1 < n2 (ramp_length=%d, n1=%d, n2=%d)",
                    traps_design.ramp_length, traps_design.n1, traps_design.n2);
  }

  /*n2 and n3 can be identical*/
  if (traps_design.n2 > traps_design.n3){
    return KS_THROW("n2 must be <= n3 (n2=%d, n3=%d)", traps_design.n2, traps_design.n3);
  }

  if (traps_design.n3 >= traps_design.n4){
    return KS_THROW("n3 must be < n4 (n3=%d, n4=%d)", traps_design.n3, traps_design.n4);
  }

  if (traps_design.n4 > etl){
    return KS_THROW("n4 must be <= etl (n4=%d, etl=%d)", traps_design.n4, etl);
  }

  if (scaling_factor < 0 || scaling_factor > 1){
    return KS_THROW("scaling factor %f out of bound [0 1]", scaling_factor);
  }

  int i;
  for (i=0; i<traps_design.ramp_length; i++){
    flip_angles[i] = (traps_design.an1+ramp[i]*(180.0-traps_design.an1))*scaling_factor;
  }

  for (i=traps_design.ramp_length; i<traps_design.an1; i++){
    flip_angles[i] = traps_design.an1*scaling_factor;
  }

  /*acending sine from an1 to an23*/
  double radi;
  int n_sine_a = traps_design.n2 - traps_design.n1;
  if (n_sine_a == 1){
    n_sine_a = 0; /* avoid devision by 0 in the loop below, in this case the value of n_sine_a is irrelevant as (i-n1) will be zero */
  }

  for (i=traps_design.n1; i<traps_design.n2; i++){
    radi = (-PI/2.0) + ((PI)/(n_sine_a-1))*(i-traps_design.n1);
    flip_angles[i] = ((sin(radi)+1)/2 * (traps_design.an23-traps_design.an1) + traps_design.an1)*scaling_factor;
  }

  /* Plateau, if n2 == n3 --> no plateau */
  for (i=traps_design.n2; i<traps_design.n3; i++){
    flip_angles[i] = traps_design.an23;
  }

  /*decending sine*/
  int n_sine_d = traps_design.n4 - traps_design.n3;
  if (n_sine_d == 1){
    n_sine_d = 0; /* avoid devision by 0 in the loop below, in this case the value of n_sine_a is irrelevant as (n3 - i) will be zero */
  }
  for (i=traps_design.n3; i<traps_design.n4; i++){
    radi = (PI/2.0) + ((PI)/(n_sine_d-1))*(traps_design.n3-i);
    flip_angles[i] = ((sin(radi)+1)/2 * (traps_design.an23-traps_design.an4) + traps_design.an4)*scaling_factor;
  }

  /* tail */
  for (i=traps_design.n4; i<etl; i++){
    flip_angles[i] = traps_design.an4*scaling_factor;
  }

  /*KS_THROW("etl=%d, ramp_length=%d, n1=%d, n2=%d, n3=%d, n4=%d, an1=%0.1f, an23=%0.1f, an4=%0.1f", etl, ramp_length, n1, n2, n3, n4, an1, an23, an4);*/
  /*
  for (i=0; i<etl; i++){
    KS_THROW("flip(%d)=%0.2f",i,flip_angles[i]);
  }
  */

  return SUCCESS; 

} /* ksepg_sineTRAPS2() */

ksepg_calc_relSAR()

double ksepg_calc_relSAR	(	double *	flip_angles,
		int	etl
	)

Computes the SAR value relative to a train with just 180° refocusing pulses

Parameters

[in]	flip_angles	Refocusing flip angles
[in]	etl	Number of echoes in the train

Returns

double Relative SAR

                                                       {

  int i; 
  double sumOfFAs = 0.0;

  for (i = 0; i < etl; i++) {
    sumOfFAs += pow(flip_angles[i]/180,2);
  }

  return (1.0/(0.25+etl)) * (0.25 + sumOfFAs);

} /* ksepg_calc_relSAR() */

Variable Documentation

rhkacq_uid

int rhkacq_uid

ks_plot_filefmt

int ks_plot_filefmt

ks_plot_kstmp

int ks_plot_kstmp

ks_psdname

char ks_psdname[256]