Commit 11e3cede authored by Cyril Poupon's avatar Cyril Poupon
Browse files

debug of SingleCompartmentRelaxometryMapperCommand; now working with multithreading

parent 148147eb
......@@ -10,7 +10,7 @@ gkg::SingleCompartmentRelaxometryMCMCFunction::
: gkg::MonteCarloMarkovChainFunction(),
_measurementCounts( measurementCounts ),
_noiseStandardDeviations( noiseStandardDeviations ),
_parameterCount( 2 )
_parameterCount( 2 + ( int32_t )measurementCounts.size() - 1 )
{
try
......
......@@ -41,6 +41,7 @@ class SingleCompartmentRelaxometryMapperCommand : public Command,
const std::vector< std::string >& stringParameters,
const std::string& fileNameOutputProtonDensityVolume,
const std::string& fileNameOutputRelaxationTimeVolume,
const std::string& fileNameOutputScalingVolume,
const std::string& fileNameFittedMeasurementVolume,
bool ascii,
const std::string& format,
......@@ -69,6 +70,7 @@ class SingleCompartmentRelaxometryMapperCommand : public Command,
const std::vector< std::string >& stringParameters,
const std::string& fileNameOutputProtonDensityVolume,
const std::string& fileNameOutputRelaxationTimeVolume,
const std::string& fileNameOutputScalingVolume,
const std::string& fileNameFittedMeasurementVolume,
bool ascii,
const std::string& format,
......
......@@ -62,7 +62,7 @@ void gkg::SingleCompartmentRelaxometryMapperGauge::display()
std::cout << " [ " << std::setw( 10 ) << _count
<< " / " << std::setw( 10 ) << _maximumCount
<< " ] " << std::flush;
<< " ]" << std::flush;
_siteIndexAlreadyDisplayed = true;
......
......@@ -28,6 +28,7 @@ gkg::SingleCompartmentRelaxometryMapperLoopContext::
const bool& verbose,
gkg::Volume< float >& protonDensityVolume,
gkg::Volume< float >& relaxationTimeVolume,
gkg::Volume< float >& scalingVolume,
const bool& computeFittedMeasurements,
gkg::Volume< float >& fittedMeasurementVolume )
: gkg::LoopContext( &gauge ),
......@@ -43,6 +44,7 @@ gkg::SingleCompartmentRelaxometryMapperLoopContext::
_verbose( verbose ),
_protonDensityVolume( protonDensityVolume ),
_relaxationTimeVolume( relaxationTimeVolume ),
_scalingVolume( scalingVolume ),
_computeFittedMeasurements( computeFittedMeasurements ),
_fittedMeasurementVolume( fittedMeasurementVolume )
{
......@@ -62,47 +64,107 @@ gkg::SingleCompartmentRelaxometryMapperLoopContext::
// initializing parameter boundaries, step and fixed information
////////////////////////////////////////////////////////////////////////////
_initialParameters.reallocate( 2 );
_initialParameters.reallocate( 2 + _inputVolumeCount - 1 );
_initialParameters( 0 ) = _scalarParameters[ 0 ];
_initialParameters( 1 ) = _scalarParameters[ 1 ];
int32_t s = 0;
for ( s = 1; s < _inputVolumeCount; s++ )
{
_initialParameters( 1 + s ) = _scalarParameters[ 1 + s ];
}
_lowerParameterBoundaries.reallocate( 2 + _inputVolumeCount - 1 );
_lowerParameterBoundaries( 0 ) = _scalarParameters[ 2 +
_inputVolumeCount - 1 ];
_lowerParameterBoundaries( 1 ) = _scalarParameters[ 3 +
_inputVolumeCount - 1 ];
for ( s = 1; s < _inputVolumeCount; s++ )
{
_lowerParameterBoundaries.reallocate( 2 );
_lowerParameterBoundaries( 0 ) = _scalarParameters[ 2 ];
_lowerParameterBoundaries( 1 ) = _scalarParameters[ 3 ];
_lowerParameterBoundaries( 1 + s ) =
_scalarParameters[ 3 + _inputVolumeCount - 1 + s ];
_upperParameterBoundaries.reallocate( 2 );
_upperParameterBoundaries( 0 ) = _scalarParameters[ 4 ];
_upperParameterBoundaries( 1 ) = _scalarParameters[ 5 ];
_deltaParameters.reallocate( 2 );
_deltaParameters( 0 ) = _scalarParameters[ 6 ];
_deltaParameters( 1 ) = _scalarParameters[ 7 ];
}
_upperParameterBoundaries.reallocate( 2 + _inputVolumeCount - 1 );
_upperParameterBoundaries( 0 ) = _scalarParameters[
4 + 2 * ( _inputVolumeCount - 1 ) ];
_upperParameterBoundaries( 1 ) = _scalarParameters[
5 + 2 * ( _inputVolumeCount - 1 ) ];
for ( s = 1; s < _inputVolumeCount; s++ )
{
_upperParameterBoundaries( 1 + s ) =
_scalarParameters[ 5 + 2 * ( _inputVolumeCount - 1 ) + s ];
_isFixedParameters.resize( 2U );
_isFixedParameters[ 0 ] = false;
_isFixedParameters[ 1 ] = false;
}
_deltaParameters.reallocate( 2 + _inputVolumeCount - 1 );
_deltaParameters( 0 ) = _scalarParameters[
6 + 3 * ( _inputVolumeCount - 1 ) ];
_deltaParameters( 1 ) = _scalarParameters[
7 + 3 * ( _inputVolumeCount - 1 ) ];
for ( s = 1; s < _inputVolumeCount; s++ )
{
_deltaParameters( 1 + s ) =
_scalarParameters[ 7 + 3 * ( _inputVolumeCount - 1 ) + s ];
}
_isFixedParameters.resize( 2U + _inputVolumeCount - 1 );
_isFixedParameters[ 0 ] = false; // force estimation of proton density
_isFixedParameters[ 1 ] = false; // force estimation of relaxation time
for ( s = 1; s < _inputVolumeCount; s++ )
{
_isFixedParameters[ 1 + s ] =
_scalarParameters[ 7 + 4 * ( _inputVolumeCount - 1 ) + s ];
}
////////////////////////////////////////////////////////////////////////////
// preparing NLP constraints and optimizer parameters
////////////////////////////////////////////////////////////////////////////
_constraints.resize( 2U );
_nlpOptimizerParameters.resize( 4U );
_constraints.resize( 2U + _inputVolumeCount - 1 );
_nlpOptimizerParameters.resize( 4U + _inputVolumeCount - 1 );
_nlpOptimizerParameters[ 0 ] = _optimizerParameters[ 0 ];
_nlpOptimizerParameters[ 1 ] = _optimizerParameters[ 1 ];
int32_t p = 0;
for ( p = 0 ; p < 2; p++ )
for ( p = 0 ; p < 2 + _inputVolumeCount - 1; p++ )
{
_constraints[ p ].reset( new gkg::ClosedIntervalOptimizerConstraint(
_lowerParameterBoundaries( p ),
_upperParameterBoundaries( p ) ) );
_nlpOptimizerParameters[ 2 + p ] =
if ( _isFixedParameters[ p ] )
{
_nlpOptimizerParameters[ 2 + p ] =
_constraints[ p ]->constrainedDeltaParameter( 0.0 );
}
else
{
_nlpOptimizerParameters[ 2 + p ] =
_constraints[ p ]->constrainedDeltaParameter( _deltaParameters( p ) );
}
}
......@@ -154,7 +216,8 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
////////////////////////////////////////////////////////////////////////////
std::vector< double >
noiseStandardDeviations( _scalarParameters.begin() + 8U,
noiseStandardDeviations( _scalarParameters.begin() +
8U + 5 * ( _inputVolumeCount - 1 ),
_scalarParameters.end() );
gkg::RCPointer< gkg::SingleCompartmentRelaxometryMCMCFunction >
......@@ -201,7 +264,21 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
std::vector< float > protonDensities( count );
std::vector< float > relaxationTimes( count );
std::vector< gkg::Vector > fittedMeasurements( count );
std::vector< std::vector< float > > scalings;
if ( _inputVolumeCount > 1 )
{
scalings.resize( count );
}
std::vector< gkg::Vector > fittedMeasurements;
if ( _computeFittedMeasurements )
{
fittedMeasurements.resize( count );
}
////////////////////////////////////////////////////////////////////////////
// looping over site(s)
......@@ -215,10 +292,11 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
double acceptanceRate = 0.0;
gkg::Matrix parameterEstimates;
gkg::MedianFilter< gkg::Vector, double > medianFilter;
gkg::Vector optimalParameters( 2 );
gkg::Vector optimalParameters( 2 + _inputVolumeCount - 1 );
std::vector< double > parameterValues( _mcmcSampleCount );
int32_t p = 0;
int32_t e = 0;
int32_t i = 0;
int32_t m = 0;
int32_t s = 0;
......@@ -326,7 +404,7 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
// computing median values for each parameter distribution
////////////////////////////////////////////////////////////////////////
for ( p = 0; p < 2; p++ )
for ( p = 0; p < 2 + _inputVolumeCount - 1; p++ )
{
for ( e = 0; e < _mcmcSampleCount; e++ )
......@@ -343,7 +421,7 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
else
{
for ( p = 0; p < 2; p++ )
for ( p = 0; p < 2 + _inputVolumeCount - 1; p++ )
{
optimalParameters( p ) = copyOfInitialParameters( p );
......@@ -359,7 +437,18 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
protonDensities[ s ] = ( float )optimalParameters( 0 );
relaxationTimes[ s ] = ( float )optimalParameters( 1 );
if ( _inputVolumeCount > 1 )
{
scalings[ s ].resize( _inputVolumeCount - 1 );
for ( i = 0; i < _inputVolumeCount - 1; i++ )
{
scalings[ s ][ i ] = ( float )optimalParameters( 2 + i );
}
}
//////////////////////////////////////////////////////////////////////////
// computing fitted data
......@@ -387,10 +476,21 @@ void gkg::SingleCompartmentRelaxometryMapperLoopContext::doIt(
for ( s = 0; s < count; s++ )
{
const gkg::Vector3d< int32_t >& site = _sites[ s ];
const gkg::Vector3d< int32_t >& site = _sites[ startIndex + s ];
_protonDensityVolume( site ) = protonDensities[ s ];
_relaxationTimeVolume( site ) = relaxationTimes[ s ];
if ( _inputVolumeCount > 1 )
{
for ( i = 0; i < _inputVolumeCount - 1; i++ )
{
_scalingVolume( site, i ) = scalings[ s ][ i ];
}
}
if ( _computeFittedMeasurements )
{
......
......@@ -40,6 +40,7 @@ class SingleCompartmentRelaxometryMapperLoopContext : public LoopContext
const bool& verbose,
Volume< float >& protonDensityVolume,
Volume< float >& relaxationTimeVolume,
Volume< float >& scalingVolume,
const bool& computeFittedMeasurements,
Volume< float >& fittedMeasurementVolume );
virtual ~SingleCompartmentRelaxometryMapperLoopContext();
......@@ -60,6 +61,7 @@ class SingleCompartmentRelaxometryMapperLoopContext : public LoopContext
const bool& _verbose;
Volume< float >& _protonDensityVolume;
Volume< float >& _relaxationTimeVolume;
Volume< float >& _scalingVolume;
const bool& _computeFittedMeasurements;
Volume< float >& _fittedMeasurementVolume;
......
......@@ -78,9 +78,6 @@ double gkg::SingleCompartmentRelaxometryNLPFunction::getLogLikelihood(
const double* currentMeasurementPointer = &_currentMeasurements( 0 );
const double* realMeasurementPointer = &_realMeasurements( 0 );
const std::vector< int32_t >& measurementCounts =
_singleCompartmentRelaxometryMCMCFunction->getMeasurementCounts();
const std::vector< double >& noiseVariances =
......
......@@ -13,7 +13,7 @@ gkg::T1IRSERelaxometryMCMCFunction::T1IRSERelaxometryMCMCFunction(
: gkg::SingleCompartmentRelaxometryMCMCFunction( measurementCounts,
noiseStandardDeviations ),
_acquisitionParameterSet(
static_cast< gkg::T1IRSEAcquisitionParameterSet* >(
static_cast< const gkg::T1IRSEAcquisitionParameterSet* >(
acquisitionParameterSet.getPointer() ) )
{
}
......@@ -41,11 +41,19 @@ void gkg::T1IRSERelaxometryMCMCFunction::getValuesAt(
int32_t m = 0;
int32_t inputVolumeIndex = 0;
int32_t measurementIndex = 0;
double scaling = 1.0;
for ( inputVolumeIndex = 0;
inputVolumeIndex < inputVolumeCount;
inputVolumeIndex++ )
{
if ( inputVolumeIndex > 0 )
{
scaling = parameters( 1 + inputVolumeIndex );
}
for ( measurementIndex = 0;
measurementIndex < _measurementCounts[ inputVolumeIndex ];
measurementIndex++ )
......@@ -53,7 +61,7 @@ void gkg::T1IRSERelaxometryMCMCFunction::getValuesAt(
values( m ) =
std::abs(
protonDensity *
scaling * protonDensity *
( 1.0 - 2 *
std::exp( -inversionTimes[ inputVolumeIndex]( measurementIndex ) /
quantitativeT1 ) ) );
......
......@@ -42,7 +42,7 @@ class T1IRSERelaxometryMCMCFunction :
const std::vector< double >&,
const RCPointer< AcquisitionParameterSet >& >;
RCPointer< T1IRSEAcquisitionParameterSet > _acquisitionParameterSet;
const T1IRSEAcquisitionParameterSet* _acquisitionParameterSet;
};
......
......@@ -13,7 +13,7 @@ gkg::T1VFASPGRRelaxometryMCMCFunction::T1VFASPGRRelaxometryMCMCFunction(
: gkg::SingleCompartmentRelaxometryMCMCFunction( measurementCounts,
noiseStandardDeviations ),
_acquisitionParameterSet(
static_cast< gkg::T1VFASPGRAcquisitionParameterSet* >(
static_cast< const gkg::T1VFASPGRAcquisitionParameterSet* >(
acquisitionParameterSet.getPointer() ) )
{
}
......@@ -56,6 +56,7 @@ void gkg::T1VFASPGRRelaxometryMCMCFunction::getValuesAt(
int32_t m = 0;
int32_t inputVolumeIndex = 0;
int32_t measurementIndex = 0;
double scaling = 1.0;
double expMinusTROverT1 = 0.0;
double correctedFlipAngle = 0.0;
for ( inputVolumeIndex = 0;
......@@ -63,6 +64,13 @@ void gkg::T1VFASPGRRelaxometryMCMCFunction::getValuesAt(
inputVolumeIndex++ )
{
if ( inputVolumeIndex > 0 )
{
scaling = parameters( 1 + inputVolumeIndex );
}
for ( measurementIndex = 0;
measurementIndex < _measurementCounts[ inputVolumeIndex ];
measurementIndex++ )
......@@ -74,7 +82,7 @@ void gkg::T1VFASPGRRelaxometryMCMCFunction::getValuesAt(
correctedFlipAngle = flipAngles[ inputVolumeIndex ]( measurementIndex )
* b1Scaling;
values( m ) = protonDensity *
values( m ) = scaling * protonDensity *
( ( 1.0 - expMinusTROverT1 ) *
std::sin( correctedFlipAngle ) ) /
( 1.0 - expMinusTROverT1 *
......
......@@ -43,7 +43,7 @@ class T1VFASPGRRelaxometryMCMCFunction :
const std::vector< double >&,
const RCPointer< AcquisitionParameterSet >& >;
RCPointer< T1VFASPGRAcquisitionParameterSet > _acquisitionParameterSet;
const T1VFASPGRAcquisitionParameterSet* _acquisitionParameterSet;
};
......
......@@ -13,7 +13,7 @@ gkg::T2MSMERelaxometryMCMCFunction::T2MSMERelaxometryMCMCFunction(
: gkg::SingleCompartmentRelaxometryMCMCFunction( measurementCounts,
noiseStandardDeviations ),
_acquisitionParameterSet(
static_cast< gkg::T2MSMEAcquisitionParameterSet* >(
static_cast< const gkg::T2MSMEAcquisitionParameterSet* >(
acquisitionParameterSet.getPointer() ) )
{
}
......@@ -41,18 +41,26 @@ void gkg::T2MSMERelaxometryMCMCFunction::getValuesAt(
int32_t m = 0;
int32_t inputVolumeIndex = 0;
int32_t measurementIndex = 0;
double scaling = 1.0;
for ( inputVolumeIndex = 0;
inputVolumeIndex < inputVolumeCount;
inputVolumeIndex++ )
{
if ( inputVolumeIndex > 0 )
{
scaling = parameters( 1 + inputVolumeIndex );
}
for ( measurementIndex = 0;
measurementIndex < _measurementCounts[ inputVolumeIndex ];
measurementIndex++ )
{
values( m ) =
protonDensity *
scaling * protonDensity *
std::exp( -echoTimes[ inputVolumeIndex ]( measurementIndex ) /
quantitativeT2 );
......
......@@ -42,7 +42,7 @@ class T2MSMERelaxometryMCMCFunction :
const std::vector< double >&,
const RCPointer< AcquisitionParameterSet >& >;
RCPointer< T2MSMEAcquisitionParameterSet > _acquisitionParameterSet;
const T2MSMEAcquisitionParameterSet* _acquisitionParameterSet;
};
......
......@@ -13,7 +13,7 @@ gkg::T2StarMGRERelaxometryMCMCFunction::T2StarMGRERelaxometryMCMCFunction(
: gkg::SingleCompartmentRelaxometryMCMCFunction( measurementCounts,
noiseStandardDeviations ),
_acquisitionParameterSet(
static_cast< gkg::T2StarMGREAcquisitionParameterSet* >(
static_cast< const gkg::T2StarMGREAcquisitionParameterSet* >(
acquisitionParameterSet.getPointer() ) )
{
}
......@@ -41,18 +41,26 @@ void gkg::T2StarMGRERelaxometryMCMCFunction::getValuesAt(
int32_t m = 0;
int32_t inputVolumeIndex = 0;
int32_t measurementIndex = 0;
double scaling = 1.0;
for ( inputVolumeIndex = 0;
inputVolumeIndex < inputVolumeCount;
inputVolumeIndex++ )
{
if ( inputVolumeIndex > 0 )
{
scaling = parameters( 1 + inputVolumeIndex );
}
for ( measurementIndex = 0;
measurementIndex < _measurementCounts[ inputVolumeIndex ];
measurementIndex++ )
{
values( m ) =
protonDensity *
scaling * protonDensity *
std::exp( -echoTimes[ inputVolumeIndex]( measurementIndex ) /
quantitativeT2Star );
......
......@@ -42,7 +42,7 @@ class T2StarMGRERelaxometryMCMCFunction :
const std::vector< double >&,
const RCPointer< AcquisitionParameterSet >& >;
RCPointer< T2StarMGREAcquisitionParameterSet > _acquisitionParameterSet;
const T2StarMGREAcquisitionParameterSet* _acquisitionParameterSet;
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment