ksdiffusion.cc File Reference

#include "ksdiffusion.h"

Functions
STATUS	ksdiff_readtensorfile (KSDIFF_SCHEME diffscheme, int ndirs)
STATUS	ksdiff_eval_scheme (KSDIFF_SCHEME logical_scheme, KSDIFF_SCHEME physical_scheme, int ndirs, ksdiff_coordinate_system coordinate_system, const KS_MAT3x3f rot_l2p)
float	ksdiff_scale_vector_to_bounds (const float *vector, const float *bounds)
float	ksdiff_vector_mag (const float *vector)
float	ksdiff_scale_vector_to_iso_bound (const float *vector, const float bound)
int	ksdiff_compare_power_mag_than_x (float *amp, double mag, float *max_power_amp_scale, double max_mag)
int	ksdiff_compare_power_x_than_mag (float *amp, double mag, float *max_power_amp_scale, double max_mag)
STATUS	ksdiff_compute_scheme_amps (KSDIFF_SCHEME diff_amp_scaling, float *max_magnitude_amp, float *RMS_amp_scale, float *max_power_amp_scale, int(*amp_cmp_f)(float *, double, float *, double), const int ndirs, const ksdiff_coordinate_system coordinate_system, const ksdiff_physical_bound_type physical_bound_type, const KS_MAT3x3f rot_l2p)
float	ksdiff_getmaxb (const float *bvalues, int n)
STATUS	ksdiff_eval_design (KSDIFF_CONTROL *diff_control, const KSDIFF_CONTROL_DESIGN *design)
s64	ksdiff_schemeloop (const KSDIFF_CONTROL *diff_control, const KSSCAN_LOOP_CONTROL *orig_loop_control, KS_DYNAMIC_STATE *dynamic, KS_CORESLICETIME coreslice(const SCAN_INFO *slice_pos, KS_DYNAMIC_STATE *dynamic))
void	ksdiff_plotloop (const KSDIFF_CONTROL *diff_control, KS_CORESLICETIME coreslice(const SCAN_INFO *slice_pos, KS_DYNAMIC_STATE *dynamic), const SCAN_INFO *scan_info)

Function Documentation

ksdiff_readtensorfile()

STATUS ksdiff_readtensorfile	(	KSDIFF_SCHEME	diffscheme,
		int	ndirs
	)

                                                                  {
    int i;
    FILE *fp;
    char tmpstr[1000];
    int lineskip = TRUE;
    char matchstr[10];
    int matchlen;

#ifdef PSD_HW
    const char *tensorfile = "/usr/g/research/tensor.dat";
    const char *tensorfile2 = "/usr/g/bin/tensor.dat";
    const char *tensorfile3 = "/opt/gehc/mr_psd_dist/root/etc/tensor.dat"; /* MR30+ */
#else
    const char *tensorfile = "./tensor.dat";
    const char *tensorfile2 = "../ksmodules/tensor.dat";
    const char *tensorfile3 = "~/tensor.dat";
#endif

    if (ndirs < 1 || ndirs > 150) {
        return KS_THROW("# diff. directions must be in range [1,150]");
    }

    /* try first 'tensorfile' */ 
    fp = fopen(tensorfile, "r" );
    if (fp == NULL) { 
      /* then try 'tensorfile2' */
      fp = fopen(tensorfile2, "r" );
      if (fp == NULL) {
        /* then try 'tensorfile3' */
        fp = fopen(tensorfile3, "r" );
        if (fp == NULL) {
          return KS_THROW("Could not open file: %s or %s or %s", tensorfile, tensorfile2, tensorfile3);
        }
      }
    }

    /* what to line content look for in the file */
    sprintf(matchstr, "%d", ndirs);
    matchlen = strlen(matchstr);

    while (lineskip) {
      fgets(tmpstr, 1000, fp);
      lineskip = strncmp(matchstr, tmpstr, matchlen);
    }

    for (i = 0; i < ndirs; i++) {
      if (fgets(tmpstr, 1000, fp) == NULL) {
        KS_THROW("Error reading file: %s for #dirs = %d, direction %d", tensorfile, ndirs, i);
      }
      sscanf(tmpstr, "%f %f %f", &diffscheme[i][XGRAD], &diffscheme[i][YGRAD], &diffscheme[i][ZGRAD]);
    }

    fclose(fp);

    return SUCCESS;

} /* ksdiff_readtensorfile() */

ksdiff_eval_scheme()

STATUS ksdiff_eval_scheme	(	KSDIFF_SCHEME	logical_scheme,
		KSDIFF_SCHEME	physical_scheme,
		int	ndirs,
		ksdiff_coordinate_system	coordinate_system,
		const KS_MAT3x3f	rot_l2p
	)

                                                    {
  int i;

  if (ndirs == 0) {
    return SUCCESS;
  }

  if (ndirs > KSDIFF_MAX_SCHEME_LENGTH && ndirs < 0) {
    return KS_THROW("Requested ndirs (%d) must be in the interval [0, %d]", ndirs, KSDIFF_MAX_SCHEME_LENGTH);
  }

  float (*diffscheme)[3] = coordinate_system == KSDIFF_LOGICAL ?
    logical_scheme : physical_scheme;

  /* clear table */
  for (i = 0; i < KSDIFF_MAX_SCHEME_LENGTH; i++) {
    diffscheme[i][XGRAD] = diffscheme[i][YGRAD] = diffscheme[i][ZGRAD] = 0.0;
  }

  /* set diffusion directions (b0 volumes are already zero, since we cleared the sequence) */
  switch (ndirs) {
    case 1: {
      diffscheme[0][XGRAD] = 0;
      diffscheme[0][YGRAD] = 0;
      diffscheme[0][ZGRAD] = 1; /* Z only */
      break;
    }
    case 2:
    case 3: {

      for (i = 0; i < 3; i++) {
        diffscheme[i][i] = 1;
      }
      break;
    }
    case 4: {

      diffscheme[0][XGRAD] = 1;
      diffscheme[0][YGRAD] = 1;
      diffscheme[0][ZGRAD] = 1;

      diffscheme[1][XGRAD] = -1;
      diffscheme[1][YGRAD] = 1;
      diffscheme[1][ZGRAD] = 1;

      diffscheme[2][XGRAD] = 1;
      diffscheme[2][YGRAD] = -1;
      diffscheme[2][ZGRAD] = 1;

      diffscheme[3][XGRAD] = -1;
      diffscheme[3][YGRAD] = -1;
      diffscheme[3][ZGRAD] = 1;

      break;
    }
    case 5:
    case 6: {

      diffscheme[0][XGRAD] = 1;
      diffscheme[0][YGRAD] = 1;
      diffscheme[0][ZGRAD] = 0;

      diffscheme[1][XGRAD] = -1;
      diffscheme[1][YGRAD] = 1;
      diffscheme[1][ZGRAD] = 0;

      diffscheme[2][XGRAD] = 1;
      diffscheme[2][YGRAD] = 0;
      diffscheme[2][ZGRAD] = 1;

      diffscheme[3][XGRAD] = -1;
      diffscheme[3][YGRAD] = 0;
      diffscheme[3][ZGRAD] = 1;

      diffscheme[4][XGRAD] = 0;
      diffscheme[4][YGRAD] = 1;
      diffscheme[4][ZGRAD] = 1;

      diffscheme[5][XGRAD] = 0;
      diffscheme[5][YGRAD] = -1;
      diffscheme[5][ZGRAD] = 1;

      break;
    }
    default: {

      STATUS s = ksdiff_readtensorfile(diffscheme, ndirs);
      KS_RAISE(s);
    }
  }

  for (i = 0; i < ndirs; i++) {
    if (coordinate_system == KSDIFF_LOGICAL) {
      ks_mat3f_apply(physical_scheme[i], rot_l2p, logical_scheme[i]);
    } else {
      ks_mat3f_invapply(logical_scheme[i], rot_l2p, physical_scheme[i]);
    }
  }

  return SUCCESS;

} /* ksdiff_eval_scheme() */

ksdiff_scale_vector_to_bounds()

float ksdiff_scale_vector_to_bounds	(	const float *	vector,
		const float *	bounds
	)

                                                         {
  float scale = 1.0f/0.0f; /* +Inf */
  int i;
  for (i=0; i<3; i++) {
    /* TODO: protect against div by zero */
    const float ratio = fabs(bounds[i]) / fabs(vector[i]);
    if (scale > ratio) {
      scale = ratio;
    }
  }
  return scale;

} /* ksdiff_scale_vector_to_bounds() */

ksdiff_vector_mag()

float ksdiff_vector_mag

(

const float *

vector

)

                                             {

  double acc = 0.0; 
  int i;
  for (i=0; i<3; i++) {
    acc += vector[i] * vector[i];
  }
  return sqrt(acc);

} /* ksdiff_vector_mag() */

ksdiff_scale_vector_to_iso_bound()

float ksdiff_scale_vector_to_iso_bound	(	const float *	vector,
		const float	bound
	)

                                                          {

  /* TODO: protect against div by zero */
  return fabs(bound) / ksdiff_vector_mag(vector);

} /* ksdiff_scale_vector_to_iso_bound() */

ksdiff_compare_power_mag_than_x()

int ksdiff_compare_power_mag_than_x	(	float *	amp,
		double	mag,
		float *	max_power_amp_scale,
		double	max_mag
	)

                                                    {

  if (mag >= max_mag &&
      amp[XGRAD] > max_power_amp_scale[XGRAD]) {
    return 1;
  }
  return 0;

} /* ksdiff_compare_power_mag_than_x() */

ksdiff_compare_power_x_than_mag()

int ksdiff_compare_power_x_than_mag	(	float *	amp,
		double	mag,
		float *	max_power_amp_scale,
		double	max_mag
	)

                                                    {

  if (amp[XGRAD] >= max_power_amp_scale[XGRAD] &&
      mag > max_mag) {
    return 1;
  }
  return 0;

} /* ksdiff_compare_power_x_than_mag() */

ksdiff_compute_scheme_amps()

STATUS ksdiff_compute_scheme_amps	(	KSDIFF_SCHEME	diff_amp_scaling,
		float *	max_magnitude_amp,
		float *	RMS_amp_scale,
		float *	max_power_amp_scale,
		int(*)(float *, double, float *, double)	amp_cmp_f,
		const int	ndirs,
		const ksdiff_coordinate_system	coordinate_system,
		const ksdiff_physical_bound_type	physical_bound_type,
		const KS_MAT3x3f	rot_l2p
	)

                                                            {
  STATUS status;
  int i, b;

  for (b = 0; b < 3; b++) {
    RMS_amp_scale[b] = 0;
    max_power_amp_scale[b] = 0;
  }

  if (ndirs < 1) {
    return SUCCESS;
  }

  /* Generate the diffusion scheme */
  /* TODO: dynamic allocation maybe? */
  KSDIFF_SCHEME logical_scheme;
  KSDIFF_SCHEME physical_scheme;
  status = ksdiff_eval_scheme(logical_scheme,
                              physical_scheme,
                              ndirs,
                              coordinate_system,
                              rot_l2p);
  KS_RAISE(status);

  /* TODO: should we check also slewrate? */
  const float logical_amp_bounds[3] = {loggrd.tx, loggrd.ty, loggrd.tz};
  const float physical_amp_bounds[3] = {phygrd.xfs, phygrd.yfs, phygrd.zfs};
  const float physical_amp_iso_bound = FMin(3, phygrd.xfs, phygrd.yfs, phygrd.zfs);

  /* Note that both scale and global_scale turn diffusion directions, whatever
     their units, into gradient amplitudes [G/cm] */
  float scale;
  float global_scale = 1.0f/0.0f; /* +Inf */
  for (i = 0; i < ndirs; i++) {

    /* Check logical constraints */
    scale = ksdiff_scale_vector_to_bounds(logical_scheme[i], logical_amp_bounds);
    if (scale < global_scale) { global_scale = scale; }

    /* Check physical constraints */
    if (physical_bound_type == KSDIFF_ISO) {      
      scale = ksdiff_scale_vector_to_iso_bound(physical_scheme[i], physical_amp_iso_bound);
    } else {
      scale = ksdiff_scale_vector_to_bounds(physical_scheme[i], physical_amp_bounds);
    }
    if (scale < global_scale) { global_scale = scale; }
  }

  /*
    now that we have the maximum feasible scaling factor we can compute the
    amp scaling factors for each direction and the absolute maximum 
    magnitude of the gradient amplitudes over all directions.
   */
  *max_magnitude_amp = 0; /* [G/cm] */

  double max_power_mag = 0; /* [G/cm] */

  if (!amp_cmp_f) {
    /* Default comparator */
    amp_cmp_f = ksdiff_compare_power_x_than_mag;
  }

  for (i = 0; i < ndirs; i++) {
    double magnitude_amp = 0;
    for (b = 0; b < 3; b++) {

      const float amp = logical_scheme[i][b]*global_scale; /* [G/cm] */
      magnitude_amp += amp * amp;

      /* belongs in [-1, 1] */
      diff_amp_scaling[i][b] = amp / logical_amp_bounds[b];

      /* Accumulate RMS amp */
      RMS_amp_scale[b] += diff_amp_scaling[i][b] * diff_amp_scaling[i][b]; 
    }
    magnitude_amp = sqrt(magnitude_amp); /* [G/cm] */
    if (magnitude_amp > *max_magnitude_amp) {
      *max_magnitude_amp = magnitude_amp;
    }

    /* Update maximum power amp */
    if (amp_cmp_f(diff_amp_scaling[i], magnitude_amp, max_power_amp_scale, max_power_mag)) {
      for (b = 0; b < 3; b++) {
        max_power_amp_scale[b] = diff_amp_scaling[i][b];
        max_power_mag = magnitude_amp;
      }
    }
  }

  /* Finalize RMS amp */
  for (b = 0; b < 3; b++) {
    RMS_amp_scale[b] = sqrt(RMS_amp_scale[b]/ndirs);
  }

  return SUCCESS;

} /* ksdiff_compute_scheme_amps() */

ksdiff_getmaxb()

float ksdiff_getmaxb	(	const float *	bvalues,
		int	n
	)

                                                  {

  int i;
  float maxb = 0.0;
  for (i = 0; i < n; i++) {
    if (bvalues[i] > maxb) {
      maxb = fabs(bvalues[i]);
    }
  }
  return maxb;

} /* ksdiff_getmaxb() */

ksdiff_eval_design()

STATUS ksdiff_eval_design	(	KSDIFF_CONTROL *	diff_control,
		const KSDIFF_CONTROL_DESIGN *	design
	)

Generate a diffusion control from its design

This function computes the scaling factors needed for the desired diffusion scheme to be executed taking into accout the limits of the gradient system according to design->physical_bound_type. The actual shape and timing of the diffusion gradients are not needed to be known, in fact, this function must be called before setting up the module that includes the diffusion gradients. In particular diff_control->bvalue and diff_control->max_magnitude_amp are to be used as inputs to design the diffusion gradients (see below for an example). Additionally, diff_control->RMS_amp_scale and diff_control->max_power_amp_scale are to be used for setting the diffusion gradients scaling factors for gradient heating calculations for average and maximum power respectively.

NOTE* that it is assumed that the diffusion gradients are generated with a maximum amplitude equal for each board equal to that board's maximum amplitude, possibly scaled down by a common factor. Each board's mximum amplitude is stored in tx, ty and tz from loggrd. In particular, one waveform (trapezoid or wave) per board is needed since those values are not necessarily equal. Specifically, the maximum amplitude of the diffusion gradients on the X logical board should be equal to alpha*loggrd.tx where alpha is a scaling factor (0 < alpha <= 1). Similarly, the maximum amplitude of the diffusion gradients on the Y and Z logical boards should be equal to alpha*loggrd.ty and alpha*loggrd.tz respectively.

For instance, in the case of standard diffusion gradients, the maximum amplitude is achieved during the plateau of the trapezoids. In order to design the shape of those trapezoids, a ramp time must be chosen first. If in doubt ks_syslimits_ramptimemax(loggrd) is always feasible as it represents the time for the slowest board to reach its maximum amplitude. Then the plateau time and the time separation between the two trapezoids should be chosen so that a plateau amplitude equal to alpha*diff_control->max_magnitude_amp results in a b-value equal to diff_control->bvalue. Finally, the actual trapezoids for the X, Y and Z boards need to be generated with a plateau amplitude equal to alpha*loggrd.tx, alpha*loggrd.ty and alpha*loggrd.tz respectively.

Parameters

[out]	diff_control
[in]	design

Returns

STATUS

                                                               {

  STATUS status;
  int i;

  /* TODO: We may want to check for Premier gradients */

  /* Copy fields */
  diff_control->nb0 = design->nb0;
  diff_control->ndirs = design->ndirs;
  diff_control->numbvals = design->numbvals;

  const float maxb = ksdiff_getmaxb(design->bvalues, design->numbvals);
  diff_control->bvalue = maxb;
  for (i=0; i < KSDIFF_MAXBVALS; ++i) {
    diff_control->multib_scale[i] = fabs(design->bvalues[i]) / maxb;
    diff_control->multib_nex[i] = design->multib_nex[i];
  }

  status = ksdiff_compute_scheme_amps(diff_control->diff_amp_scaling,
                                      &diff_control->max_magnitude_amp,
                                      diff_control->RMS_amp_scale,
                                      diff_control->max_power_amp_scale,
                                      NULL, /* Default comparator */
                                      design->ndirs,
                                      design->coordinate_system,
                                      design->physical_bound_type,
                                      design->rot_l2p);
  KS_RAISE(status);

  return SUCCESS;

} /* ksdiff_eval_design() */

ksdiff_schemeloop()

s64 ksdiff_schemeloop	(	const KSDIFF_CONTROL *	diff_control,
		const KSSCAN_LOOP_CONTROL *	loop_control,
		KS_DYNAMIC_STATE *	dynamic,
		KS_CORESLICETIME	coresliceconst SCAN_INFO *slice_pos, KS_DYNAMIC_STATE *dynamic
	)

Executes all loops for a diffusion scheme according to the provided loop control

The index dynamic->vol is incremented for each b0 and each diffusion volumes. This corresponds to a total number of sampled volumes equal to nb0 + numbvals*ndirs for each invocation of ksdiff_schemeloop.

Note that dynamic->vol is not reset to zero so that, if the diffusion scheme is repeated, all volumes have different index.

Parameters

[in]	diff_control
[in]	loop_control
[in,out]	dynamic
[in]	coreslice

Returns

s64 Time necessary for performing all relevant loops on the hardware

                                                                                                         {

  s64 time = 0;

  KS_DYNAMIC_STATE default_dynamic = KS_INIT_DYNAMIC_STATE;
  if (!dynamic) {
    dynamic = &default_dynamic;
  }  

  /* make a copy of the loop control as we are going to modify some */
  KSSCAN_LOOP_CONTROL loop_control = *orig_loop_control;

  const int numbval_indices = diff_control->nb0 + diff_control->numbvals;
  int board;

  /* Initialize for b0s */
  loop_control.naverages = orig_loop_control->naverages;
  float bvalue_scale = 0;
  int ndirs = 1;

  for (dynamic->b_idx = 0; dynamic->b_idx < numbval_indices; dynamic->b_idx++) {

    if (dynamic->b_idx >= diff_control->nb0) {
      const int bvalue_idx = dynamic->b_idx - diff_control->nb0;
      loop_control.naverages = diff_control->multib_nex[bvalue_idx];
      bvalue_scale = diff_control->multib_scale[bvalue_idx];
      ndirs = diff_control->ndirs;
    }

    for (dynamic->diff_idx = 0; dynamic->diff_idx < ndirs; dynamic->diff_idx++) {
      for (dynamic->pass = 0; dynamic->pass < loop_control.slicetiming.slice_plan.npasses; dynamic->pass++) {

        /* In a single pass scan we should play dummies only once */
        const int play_dda = loop_control.slicetiming.slice_plan.npasses == 1 ?
          (!(dynamic->vol) && !dynamic->b_idx && !dynamic->diff_idx) : 1;
        loop_control.dda = play_dda * orig_loop_control->dda;

        /*
          Set up the diffusion amplitudes

          Note: diff_amp assumes that the diffusion gradients are generated
          with amplitude equal to loggrd.t*
        */
        for (board = 0; board < 3; board++) {
          dynamic->diff_amp[board] = bvalue_scale > 0 ?
            bvalue_scale*diff_control->diff_amp_scaling[dynamic->diff_idx][board] : 0.0f;
        }

        time += ksscan_acqloop(&loop_control, dynamic, coreslice);

        loop_control.dda = 0;
        if (loop_control.play_endofpass) {
          time += GEReq_endofpass();
        }
        ks_plot_slicetime_endofpass(KS_PLOT_PASS_WAS_STANDARD);

      } /* pass loop */

      /* Done with a volume */
      dynamic->vol++;

    } /* diff_idx loop */
  } /* b_idx loop */

  return time; /* in [us] */

} /* ksdiff_schemeloop() */

ksdiff_plotloop()

void ksdiff_plotloop	(	const KSDIFF_CONTROL *	diff_control,
		KS_CORESLICETIME	coresliceconst SCAN_INFO *slice_pos, KS_DYNAMIC_STATE *dynamic,
		const SCAN_INFO *	scan_info
	)

Executes a loop over all diffusion direction for plotting purposes

All loops other than the one over all diffusion directions are collapsed.

Parameters

[in]	diff_control
[in]	coreslice
[in]	scan_info	Optional slice position, can be NULL

Returns

void

                                                 {

  KS_DYNAMIC_STATE dynamic = KS_INIT_DYNAMIC_STATE;

  SCAN_INFO default_scan_info = DEFAULT_AXIAL_SCAN_INFO;
  if (!scan_info) {
    scan_info = &default_scan_info;
  }

  /* Once for b0 */
  coreslice(scan_info, &dynamic);

  int dir, board;
  for (dir = 0; dir < diff_control->ndirs; dir++) {
    for (board = 0; board < 3; board++) {
      dynamic.diff_amp[board] =
        diff_control->diff_amp_scaling[dir][board];
    }

    coreslice(scan_info, &dynamic);
  }

} /* ksdiff_plotloop() */