Commit 148147eb authored by Cyril Poupon's avatar Cyril Poupon
Browse files

added a value to compute PGSE little & big delta from gradient characteristics and b-value target

parent 64c943e3
......@@ -456,6 +456,7 @@ gkg-dmri-qspace-sampling/SmsVb15MultipleDifferentShellQSpaceSampling.h
gkg-dmri-qspace-sampling/SmsVb15MultipleSameShellQSpaceSampling.h
gkg-dmri-qspace-sampling/SmsVb15SingleShellQSpaceSampling.h
gkg-dmri-qspace-sampling/SphericalQSpaceSampling.h
gkg-dmri-sequence-tuning/PGSESequenceTuning.h
gkg-dmri-sh-basis/SymmetricalSphericalHarmonicsCache.h
gkg-dmri-simulator-membrane/Membrane.h
gkg-dmri-simulator-motion/MotionModelFactory.h
......@@ -883,6 +884,7 @@ gkg-dmri-qspace-sampling/SmsVb15MultipleDifferentShellQSpaceSampling.cxx
gkg-dmri-qspace-sampling/SmsVb15MultipleSameShellQSpaceSampling.cxx
gkg-dmri-qspace-sampling/SmsVb15SingleShellQSpaceSampling.cxx
gkg-dmri-qspace-sampling/SphericalQSpaceSampling.cxx
gkg-dmri-sequence-tuning/PGSESequenceTuning.cxx
gkg-dmri-sh-basis/SymmetricalSphericalHarmonicsCache.cxx
gkg-dmri-simulator-membrane/Membrane.cxx
gkg-dmri-simulator-motion/MotionModel.cxx
......
#include <gkg-dmri-sequence-tuning/PGSESequenceTuning.h>
#include <gkg-core-exception/Exception.h>
#include <cmath>
#define MT_PER_M_TO_G_PER_CM 0.1
#define T_PER_M_PER_S_TO_G_PER_CM_PER_US 0.0001
#define MS_TO_US 1000.0
#define TWOPI_GAM 26748.0
#define GAM 26748.0 / ( 2 * M_PI )
#define RUP_GRD( A, GRAD_UPDATE_TIME ) \
( ( GRAD_UPDATE_TIME - ( ( int32_t )( A ) % \
GRAD_UPDATE_TIME ) ) % GRAD_UPDATE_TIME + \
( int32_t )( A ) )
/*
//
// Stejskal-Tanner' equation for trapezoidal gradients:
//
// targetAmp ----> +--------------+ +--------------+
// / \ / \
// / \ / \
// ---+--+--------------+--+------+--+--------------+--+--------
//
// rampWidth ---<->-
// pulseWidth timeSeparation
// <--------------------><----->
//
*/
void gkg::getPGSEParameters( double maximumGradientAmplitude, // in mT/m
double minimumSlewRate, // in T/m/s
double gradientResolution, // in us
double timeSeparation, // in ms
double timeBeforeEcho, // in ms
double bValue, // in s/mm2
double& littleDelta, // in ms
double& bigDelta, // in ms
double& rampWidth, // in ms
double& diffusionTime, // in ms
double& echoTime ) // in ms
{
try
{
maximumGradientAmplitude *= MT_PER_M_TO_G_PER_CM;
minimumSlewRate *= T_PER_M_PER_S_TO_G_PER_CM_PER_US;
timeSeparation *= MS_TO_US;
rampWidth = ( double )RUP_GRD( ( int32_t )( maximumGradientAmplitude /
minimumSlewRate ),
( int32_t )gradientResolution );
// solving the equation b = (GAM)^2*g^2*(d^2*(D - d/3.0)+r^3/30 - d*r^2/6)
// for "d" which results in a cubic equation (the r^ are the ramp times to
// account for the trapezoid diffusion lobes):
//
// d^3 - 3D*d^2 + dr^2/2B -(r^3/(10*B) -3b/beta) = 0
// where beta = (GAM*g)^2 and D = (pw_sep+d). Thus, the equ. becomes..
//
// d^3 + (3/2*pw_sep)*d^2 - (r^2/4)*d - (3*b/2beta - r^3/20) = 0
// casting 180/crusher time in seconds
double sep = ( double )( timeSeparation / 1.0e6 );
// diffusion ramps in seconds
double ramp = ( double )( rampWidth / 1.0e6 );
double beta = TWOPI_GAM * TWOPI_GAM * ( maximumGradientAmplitude / 10.0 ) *
( maximumGradientAmplitude / 10.0 );
double a = ( 3.0 / 2.0 ) * sep;
double b = ramp * ramp / 4;
double c = ( ( -3.0 / 2.0 ) * bValue ) / beta +
( ( ramp * ramp * ramp ) / 20.0 );
double Q = ( a * a - 3.0 * b ) / 9.0;
double R = ( ( 2.0 * ( a * a * a ) - 9.0 * a * b + 27.0 * c ) / 54.0 );
double Rsq = R * R;
double Qcub = Q * Q * Q;
double theta = 0.0;
double root[ 3 ];
double pulseWidth = 0.0;
if ( Rsq < Qcub )
{
theta = std::acos( R / std::sqrt( Qcub ) );
// casting roots in us
root[ 0 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
std::cos( theta / 3.0 ) - a / 3.0 );
root[ 1 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
std::cos( ( theta + 2.0 * M_PI ) / 3.0 ) - a / 3.0 );
root[ 2 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
std::cos( ( theta - 2.0 * M_PI ) / 3.0 ) - a / 3.0 );
// order roots and take the smallest - non negative value
int32_t i = 0;
int32_t j = 0;
double temp = 0.0;
for ( i = 0; i <= 2; i++ )
{
temp = root[ i ];
for ( j = i + 1; j <= 2; j++ )
{
if ( root[ j ] < temp )
{
root[ i ] = root[ j ];
root[ j ] = temp;
temp = root[ i ];
}
}
}
// take smallest that is larger than gradientResolution
if ( root[ 0 ] >= gradientResolution )
{
pulseWidth = ( double )( ( int32_t )root[ 0 ] );
}
else if ( root[ 1 ] >= gradientResolution )
{
pulseWidth = ( double )( ( int32_t )root[ 1 ] );
}
else if ( root[ 2 ] >= gradientResolution )
{
pulseWidth = ( double )( ( int32_t )root[ 2 ] );
}
else
{
throw std::runtime_error(
"pulse width would be lower than gradient resolution" );
}
}
else
{
double sign = -1.0;
// checking sign of R
if ( R >= 0 )
{
sign = 1.0;
}
double A = -sign * std::cbrt( std::fabs( R ) + std::sqrt( Rsq - Qcub ) );
double B = Q / A;
// this is the only real root - the other two are complex
root[ 0 ] = 1e6 * ( A + B - a / 3.0 );
if ( root[ 0 ] >= gradientResolution )
{
pulseWidth = ( double )( ( int32_t )root[ 0 ] );
}
else
{
throw std::runtime_error(
"pulse width would be lower than gradient resolution" );
}
}
// put pulseWidth on gradient resolution boundaries
pulseWidth = ( double )RUP_GRD( pulseWidth, ( int32_t )gradientResolution );
littleDelta = ( pulseWidth + 2 * rampWidth ) / MS_TO_US;
bigDelta = ( pulseWidth + 2 * rampWidth + timeSeparation ) / MS_TO_US;
diffusionTime = ( bigDelta - littleDelta / 3.0 );
echoTime = 2.0 * ( 0.5 * timeSeparation / MS_TO_US +
pulseWidth / MS_TO_US + timeBeforeEcho );
}
GKG_CATCH( "void gkg::getPGSEParameters( "
"double maximumGradientAmplitude, "
"double minimumSlewRate, "
"double gradientResolution, "
"double timeSeparation, "
"double timeBeforeEcho, "
"double bValue, "
"double& littleDelta, "
"double& bigDelta, "
"double& rampWidth, "
"double& diffusionTime, "
"double& echoTime )" );
}
#undef MT_PER_M_TO_G_PER_CM
#undef T_PER_M_PER_S_TO_G_PER_CM_PER_US
#undef MS_TO_US
#undef TWOPI_GAM
#undef GAM
#undef RUP_GRD
#ifndef _gkg_dmri_sequence_tuning_PGSESequenceTuning_h_
#define _gkg_dmri_sequence_tuning_PGSESequenceTuning_h_
namespace gkg
{
void getPGSEParameters( double maximumGradientAmplitude, // in mT/m
double minimumSlewRate, // in T/m/s
double gradientResolution, // in us
double timeSeparation, // in ms
double timeBeforeEcho, // in ms
double bValue, // in s/mm2
double& littleDelta, // in ms
double& bigDelta, // in ms
double& rampWidth, // in ms
double& diffusionTime, // in ms
double& echoTime ); // in ms
}
#endif
#include <gkg-dmri-plugin-functors/DwiBValue/DwiBValueCommand.h>
#include <gkg-communication-command/CommandFactory.h>
#include <gkg-communication-getopt/Application.h>
#include <gkg-dmri-sequence-tuning/PGSESequenceTuning.h>
#include <gkg-core-exception/Exception.h>
#include <iostream>
#include <cmath>
#include <iostream>
#define TWOPI_GAM 26748.0
#define GAM 26748.0/ (2*M_PI)
#define GRAD_UPDATE_TIME 4
#define RUP_GRD( A ) ( ( GRAD_UPDATE_TIME - ( ( int32_t )( A ) % \
GRAD_UPDATE_TIME ) ) % GRAD_UPDATE_TIME + \
( int32_t )( A ) )
#define MT_PER_M_TO_G_PER_CM 0.1
#define T_PER_M_PER_S_TO_G_PER_CM_PER_US 0.0001
#define MS_TO_US 1000.0
/*
Stejskal-Tanner' equation for trapezoidal gradients:
targetAmp ----> +--------------+ +--------------+
/ \ / \
/ \ / \
---+--+--------------+--+------+--+--------------+--+--------
rampWidth ---<->-
pulseWidth separation
<--------------------><----->
*/
double GkgGetDiffusionPulseWidth( double maxAmp,
double minSlewRate,
double separation,
double bValue,
double& rampWidth )
{
int32_t i, j;
double Q, R, Rsq, Qcub, theta;
double sep, ramp;
float temp;
double a, b, c;
double beta;
float root[ 3 ];
double pulseWidth = 0;
rampWidth = ( double )RUP_GRD( ( int32_t )( maxAmp / minSlewRate ) );
// solving the equation b = (GAM)^2*g^2*(d^2*(D - d/3.0)+r^3/30 - d*r^2/6)
// for "d" which results in a cubic equation (the r^ are the ramp times to
// account for the trapezoid diffusion lobes):
//
// d^3 - 3D*d^2 + dr^2/2B -(r^3/(10*B) -3b/beta) = 0
// where beta = (GAM*g)^2 and D = (pw_sep+d). Thus, the equ. becomes..
//
// d^3 + (3/2*pw_sep)*d^2 - (r^2/4)*d - (3*b/2beta - r^3/20) = 0
sep = ( double )( separation / 1.0e6 ); // cast 180/crusher time in seconds
ramp = ( double )( rampWidth / 1.0e6 ); // diffusion ramps in seconds
beta = TWOPI_GAM * TWOPI_GAM * ( maxAmp / 10.0 ) * ( maxAmp / 10.0 );
a = ( 3.0 / 2 ) * sep;
b = ramp * ramp / 4;
c = ( ( -3.0 / 2 ) * bValue ) / beta + ( ( ramp * ramp * ramp ) / 20.0 );
Q = ( double )( a * a - 3.0 * b ) / 9.0;
R = ( double )( ( 2.0 * ( a * a * a ) - 9.0 * a * b + 27 * c ) / 54.0 );
Rsq = ( double )( R * R );
Qcub = ( double )( Q * Q * Q );
if ( Rsq < Qcub )
{
theta = acos( R / sqrt( Qcub ) );
// cast roots in terms of us
root[ 0 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
cos( theta / 3.0 ) - a / 3.0 );
root[ 1 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
cos( ( theta + 2.0 * M_PI ) / 3.0 ) -
a / 3.0 );
root[ 2 ] = 1e6 * ( -2.0 * std::sqrt( Q ) *
cos( ( theta - 2.0 * M_PI ) / 3.0 ) -
a / 3.0 );
// order roots and take the smallest - non negative value
for ( i = 0; i <= 2; i++ )
{
temp = root[ i ];
for ( j = i + 1; j <= 2; j++ )
{
if ( root[ j ] < temp )
{
root[ i ] = root[ j ];
root[ j ] = temp;
temp = root[ i ];
}
}
}
// take smallest that is larger than GRAD_UPDATE_TIME
if ( root[ 0 ] > GRAD_UPDATE_TIME )
{
pulseWidth = ( double )( ( int32_t )root[ 0 ] );
}
else if ( root[ 1 ] > GRAD_UPDATE_TIME )
{
pulseWidth = ( double )( ( int32_t )root[ 1 ] );
}
else if ( root[ 2 ] > GRAD_UPDATE_TIME )
{
pulseWidth = ( double )( ( int32_t )root[ 2 ] );
}
}
else
{
double sign;
double A, B;
// check sign of R
if ( R >= 0 )
{
sign = 1.0;
}
else
{
sign = -1.0;
}
A = -sign * cbrt( std::fabs( R ) + std::sqrt( Rsq - Qcub ) );
B = Q / A;
// this is the only real root - the other two are complex
root[ 0 ] = 1e6 * ( float )( A + B - a / 3.0 );
if ( root[ 0 ] >= GRAD_UPDATE_TIME )
{
pulseWidth = ( double )( ( int32_t )root[ 0 ] );
}
}
// put pulseWidth on GRAD_UPDATE boundaries
pulseWidth = ( double )RUP_GRD( pulseWidth );
return pulseWidth;
}
#define GAM 26748.0/ ( 2 * M_PI )
//
......@@ -201,10 +37,10 @@ gkg::DwiBValueCommand::DwiBValueCommand( int32_t argc,
gkg::DwiBValueCommand::DwiBValueCommand( double maximuGradientAmplitude,
double minimumSlewRate,
double gradientResolution,
double timeSeparation,
double timeBeforeEcho,
double bValue,
double echoSpacing,
int32_t beforeEchoKSpaceLineCount,
double B0,
double ADC,
double averageT2 )
......@@ -216,10 +52,10 @@ gkg::DwiBValueCommand::DwiBValueCommand( double maximuGradientAmplitude,
execute( maximuGradientAmplitude,
minimumSlewRate,
gradientResolution,
timeSeparation,
timeBeforeEcho,
bValue,
echoSpacing,
beforeEchoKSpaceLineCount,
B0,
ADC,
averageT2 );
......@@ -228,10 +64,10 @@ gkg::DwiBValueCommand::DwiBValueCommand( double maximuGradientAmplitude,
GKG_CATCH( "gkg::DwiBValueCommand::DwiBValueCommand( "
"double maximuGradientAmplitude, "
"double minimumSlewRate, "
"double gradientResolution, "
"double timeSeparation, "
"double timeBeforeEcho, "
"double bValue, "
"double echoSpacing, "
"int32_t beforeEchoKSpaceLineCount, "
"double B0, "
"double ADC, "
"double averageT2 )" );
......@@ -248,20 +84,20 @@ gkg::DwiBValueCommand::DwiBValueCommand( const gkg::Dictionary& parameters )
DECLARE_FLOATING_PARAMETER( parameters, double, maximuGradientAmplitude );
DECLARE_FLOATING_PARAMETER( parameters, double, minimumSlewRate );
DECLARE_FLOATING_PARAMETER( parameters, double, gradientResolution );
DECLARE_FLOATING_PARAMETER( parameters, double, timeSeparation );
DECLARE_FLOATING_PARAMETER( parameters, double, timeBeforeEcho );
DECLARE_FLOATING_PARAMETER( parameters, double, bValue );
DECLARE_FLOATING_PARAMETER( parameters, double, echoSpacing );
DECLARE_INTEGER_PARAMETER( parameters, int32_t, beforeEchoKSpaceLineCount );
DECLARE_FLOATING_PARAMETER( parameters, double, B0 );
DECLARE_FLOATING_PARAMETER( parameters, double, ADC );
DECLARE_FLOATING_PARAMETER( parameters, double, averageT2 );
execute( maximuGradientAmplitude,
minimumSlewRate,
gradientResolution,
timeSeparation,
timeBeforeEcho,
bValue,
echoSpacing,
beforeEchoKSpaceLineCount,
B0,
ADC,
averageT2 );
......@@ -304,10 +140,10 @@ void gkg::DwiBValueCommand::parse()
double maximuGradientAmplitude = 0.0;
double minimumSlewRate = 0.0;
double gradientResolution = 1.0; // us
double timeSeparation = 0.0;
double timeBeforeEcho = 0;
double bValue = 0.0;
double echoSpacing = 0.0;
int32_t beforeEchoKSpaceLineCount = 0;
double B0 = 0.0;
double ADC = 0.0;
double averageT2 = 0.0;
......@@ -321,21 +157,22 @@ void gkg::DwiBValueCommand::parse()
application.addSingleOption( "-minimumSlewRate",
"Minimum slew rate in T/m/s",
minimumSlewRate );
application.addSingleOption( "-gradientResolution",
"Gradient resolution in us (default=1us)",
gradientResolution,
true );
application.addSingleOption( "-timeSeparation",
"Time separation between the two gradients "
"in ms",
timeSeparation );
application.addSingleOption( "-timeBeforeEcho",
"Time after end of second diffusion pulse and "
" before echo time in ms",
timeBeforeEcho,
true );
application.addSingleOption( "-bValue",
"b-balue in s/mm2",
bValue );
application.addSingleOption( "-echoSpacing",
"Echo spacing in ms",
echoSpacing,
true );
application.addSingleOption( "-beforeEchoKSpaceLineCount",
"Before echo k-space line count",
beforeEchoKSpaceLineCount,
true );
application.addSingleOption( "-B0",
"Static magnetic field in T",
B0,
......@@ -353,10 +190,10 @@ void gkg::DwiBValueCommand::parse()
execute( maximuGradientAmplitude,
minimumSlewRate,
gradientResolution,
timeSeparation,
timeBeforeEcho,
bValue,
echoSpacing,
beforeEchoKSpaceLineCount,
B0,
ADC,
averageT2 );
......@@ -369,10 +206,10 @@ void gkg::DwiBValueCommand::parse()
void gkg::DwiBValueCommand::execute( double maximuGradientAmplitude,
double minimumSlewRate,
double gradientResolution,
double timeSeparation,
double timeBeforeEcho,
double bValue,
double echoSpacing,
int32_t beforeEchoKSpaceLineCount,
double B0,
double ADC,
double averageT2 )
......@@ -381,73 +218,64 @@ void gkg::DwiBValueCommand::execute( double maximuGradientAmplitude,
try
{
maximuGradientAmplitude *= MT_PER_M_TO_G_PER_CM;
minimumSlewRate *= T_PER_M_PER_S_TO_G_PER_CM_PER_US;
timeSeparation *= MS_TO_US;
double littleDelta = 0.0;
double bigDelta = 0.0;
double rampWidth = 0.0;
double diffusionTime = 0.0;
double echoTime = 0.0;
gkg::getPGSEParameters( maximuGradientAmplitude,
minimumSlewRate,
gradientResolution,
timeSeparation,
timeBeforeEcho,
bValue,
littleDelta,