Commit e4b39cbf authored by Cyril Poupon's avatar Cyril Poupon
Browse files

new single compartment relaxometry mapper

parent cafda8bd
......@@ -62,6 +62,22 @@ class MedianFilter< std::vector< std::complex< IN > >, std::complex< OUT > > :
};
//
// class MedianFilter< gkg::Vector, class OUT >
//
template < class OUT >
class MedianFilter< gkg::Vector, OUT > :
public FilterFunction< gkg::Vector, OUT >
{
public:
void filter( const gkg::Vector& in, OUT& out ) const;
};
//
// class MedianFilter< Volume< IN >, OUT >
//
......
......@@ -119,6 +119,53 @@ void gkg::MedianFilter< std::vector< std::complex< IN > >,
}
//
// class MedianFilter< gkg::Vector, OUT >
//
template < class OUT >
inline
void gkg::MedianFilter< gkg::Vector, OUT >::filter(
const gkg::Vector& in, OUT& out ) const
{
try
{
if ( in.getSize() == 0 )
{
throw std::runtime_error( "empty vector" );
}
std::vector< double > copy( in.getSize() );
int32_t index = 0;
std::vector< double >::iterator
c = copy.begin(),
ce = copy.end();
while ( c != ce )
{
*c = in( index );
++ index;
++ c;
}
std::sort( copy.begin(), copy.end() );
out = ( OUT )copy[ copy.size() / 2 - ( copy.size() % 2 ? 0U : 1U ) ];
}
GKG_CATCH( "template < class OUT > "
"inline "
"void gkg::MedianFilter< gkg::Vector, OUT >::filter( "
"const gkg::Vector& in, OUT& out ) const" );
}
//
// class MedianFilter< Volume< IN >, OUT >
//
......
......@@ -20,6 +20,23 @@ MyelinWaterFractionMapper/VFASPGRBasedMyelinWaterFractionMCMCFunction.h
RelaxometryKalmanFilter/RelaxometryKalmanFilterCommand.h
RelaxometryMapper/RelaxometryMapperCommand.h
RelaxometryT1CorticalLaminationMapper/RelaxometryT1CorticalLaminationMapperCommand.h
SingleCompartmentRelaxometryMapper/AcquisitionParameterSet.h
SingleCompartmentRelaxometryMapper/AcquisitionParameterSetFactory.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperCommand.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperGauge.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperLoopContext.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunctionFactory.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunction.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCParameterSet.h
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryNLPFunction.h
SingleCompartmentRelaxometryMapper/T1IRSEAcquisitionParameterSet.h
SingleCompartmentRelaxometryMapper/T1IRSERelaxometryMCMCFunction.h
SingleCompartmentRelaxometryMapper/T1VFASPGRAcquisitionParameterSet.h
SingleCompartmentRelaxometryMapper/T1VFASPGRRelaxometryMCMCFunction.h
SingleCompartmentRelaxometryMapper/T2MSMEAcquisitionParameterSet.h
SingleCompartmentRelaxometryMapper/T2MSMERelaxometryMCMCFunction.h
SingleCompartmentRelaxometryMapper/T2StarMGREAcquisitionParameterSet.h
SingleCompartmentRelaxometryMapper/T2StarMGRERelaxometryMCMCFunction.h
)
......@@ -45,6 +62,23 @@ MyelinWaterFractionMapper/VFASPGRBasedMyelinWaterFractionMCMCFunction.cxx
RelaxometryKalmanFilter/RelaxometryKalmanFilterCommand.cxx
RelaxometryMapper/RelaxometryMapperCommand.cxx
RelaxometryT1CorticalLaminationMapper/RelaxometryT1CorticalLaminationMapperCommand.cxx
SingleCompartmentRelaxometryMapper/AcquisitionParameterSet.cxx
SingleCompartmentRelaxometryMapper/AcquisitionParameterSetFactory.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperCommand.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperGauge.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperLoopContext.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunction.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunctionFactory.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCParameterSet.cxx
SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryNLPFunction.cxx
SingleCompartmentRelaxometryMapper/T1IRSEAcquisitionParameterSet.cxx
SingleCompartmentRelaxometryMapper/T1IRSERelaxometryMCMCFunction.cxx
SingleCompartmentRelaxometryMapper/T1VFASPGRAcquisitionParameterSet.cxx
SingleCompartmentRelaxometryMapper/T1VFASPGRRelaxometryMCMCFunction.cxx
SingleCompartmentRelaxometryMapper/T2MSMEAcquisitionParameterSet.cxx
SingleCompartmentRelaxometryMapper/T2MSMERelaxometryMCMCFunction.cxx
SingleCompartmentRelaxometryMapper/T2StarMGREAcquisitionParameterSet.cxx
SingleCompartmentRelaxometryMapper/T2StarMGRERelaxometryMCMCFunction.cxx
)
......
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunction.h>
#include <gkg-core-exception/Exception.h>
#include <cmath>
gkg::SingleCompartmentRelaxometryMCMCFunction::
SingleCompartmentRelaxometryMCMCFunction(
const std::vector< int32_t >& measurementCounts,
const std::vector< double >& noiseStandardDeviations )
: gkg::MonteCarloMarkovChainFunction(),
_measurementCounts( measurementCounts ),
_noiseStandardDeviations( noiseStandardDeviations ),
_parameterCount( 2 )
{
try
{
int32_t inputVolumeCount = ( int32_t )_measurementCounts.size();
if ( ( int32_t )_noiseStandardDeviations.size() != inputVolumeCount )
{
throw std::runtime_error(
"noise standard deviation count does not match input volume count" );
}
_noiseVariances.resize( inputVolumeCount );
int32_t inputVolumeIndex = 0;
for ( inputVolumeIndex = 0;
inputVolumeIndex < inputVolumeCount;
inputVolumeIndex++ )
{
_noiseVariances[ inputVolumeIndex ] =
_noiseStandardDeviations[ inputVolumeIndex ] *
_noiseStandardDeviations[ inputVolumeIndex ];
}
}
GKG_CATCH( "gkg::SingleCompartmentRelaxometryMCMCFunction:: "
"SingleCompartmentRelaxometryMCMCFunction( "
"const std::vector< int32_t >& measurementCounts, "
"const std::double< double >& noiseStandardDeviations )" );
}
gkg::SingleCompartmentRelaxometryMCMCFunction::
~SingleCompartmentRelaxometryMCMCFunction()
{
}
double gkg::SingleCompartmentRelaxometryMCMCFunction::getLikelihoodRatio(
const gkg::Vector& realMeasurements,
const gkg::Vector& currentMeasurements,
const gkg::Vector& newMeasurements ) const
{
try
{
int32_t inputVolumeCount = ( int32_t )_measurementCounts.size();
// Gaussian case
double logLikelihood1 = 0.0;
double logLikelihood2 = 0.0;
int32_t m = 0;
int32_t measurementIndex = 0;
int32_t inputVolumeIndex = 0;
for ( inputVolumeIndex = 0;
inputVolumeIndex < inputVolumeCount;
inputVolumeIndex++ )
{
for ( measurementIndex = 0;
measurementIndex < _measurementCounts[ inputVolumeIndex ];
measurementIndex++ )
{
logLikelihood1 += -std::log( _noiseStandardDeviations[
inputVolumeIndex ] *
std::sqrt( 2 * M_PI ) ) -
( ( currentMeasurements( m ) -
realMeasurements( m ) ) *
( currentMeasurements( m ) -
realMeasurements( m ) ) ) /
( 2.0 * _noiseVariances[ inputVolumeIndex ] );
logLikelihood2 += -std::log( _noiseStandardDeviations[
inputVolumeIndex ] *
std::sqrt( 2 * M_PI ) ) -
( ( newMeasurements( m ) -
realMeasurements( m ) ) *
( newMeasurements( m ) -
realMeasurements( m ) ) ) /
( 2.0 * _noiseVariances[ inputVolumeIndex ] );
++ m;
}
}
return std::exp( logLikelihood2 - logLikelihood1 );
}
GKG_CATCH( "double gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getLikelihoodRation( "
"const gkg::Vector& realMeasurements, "
"const gkg::Vector& currentMeasurements, "
"const gkg::Vector& newMeasurements ) const" );
}
const std::vector< int32_t >&
gkg::SingleCompartmentRelaxometryMCMCFunction::getMeasurementCounts() const
{
try
{
return _measurementCounts;
}
GKG_CATCH( "const std::vector< int32_t >& "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getMeasurementCounts() const" );
}
int32_t
gkg::SingleCompartmentRelaxometryMCMCFunction::getMeasurementCount(
int32_t inputVolumeIndex ) const
{
try
{
return _measurementCounts[ inputVolumeIndex ];
}
GKG_CATCH( "int32_t "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getMeasurementCount( "
"int32_t inputVolumeIndex ) const" );
}
const std::vector< double >&
gkg::SingleCompartmentRelaxometryMCMCFunction::
getNoiseStandardDeviations() const
{
try
{
return _noiseStandardDeviations;
}
GKG_CATCH( "const std::vector< double >& "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getNoiseStandardDeviations() const" );
}
double
gkg::SingleCompartmentRelaxometryMCMCFunction::
getNoiseStandardDeviation( int32_t inputVolumeIndex ) const
{
try
{
return _noiseStandardDeviations[ inputVolumeIndex ];
}
GKG_CATCH( "double "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getNoiseStandardDeviation( "
"int32_t inputVolumeIndex ) const" );
}
const std::vector< double >&
gkg::SingleCompartmentRelaxometryMCMCFunction::getNoiseVariances() const
{
try
{
return _noiseVariances;
}
GKG_CATCH( "const std::vector< double >& "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getNoiseVariancess() const" );
}
double
gkg::SingleCompartmentRelaxometryMCMCFunction::
getNoiseVariance( int32_t inputVolumeIndex ) const
{
try
{
return _noiseVariances[ inputVolumeIndex ];
}
GKG_CATCH( "double "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getNoiseVariance( "
"int32_t inputVolumeIndex ) const" );
}
int32_t
gkg::SingleCompartmentRelaxometryMCMCFunction::getParameterCount() const
{
try
{
return _parameterCount;
}
GKG_CATCH( "int32_t "
"gkg::SingleCompartmentRelaxometryMCMCFunction::"
"getParameterCount() const" );
}
void gkg::SingleCompartmentRelaxometryMCMCFunction::setSite(
const gkg::Vector3d< int32_t >& site )
{
try
{
_site = site;
}
GKG_CATCH( "void gkg::SingleCompartmentRelaxometryMCMCFunction::setSite( "
"const gkg::Vector3d< int32_t >& site )" );
}
#ifndef _gkg_qmri_plugin_functors_SingleCompartmentRelaxometryMapper_SingleCompartmentRelaxometryMCMCFunction_h_
#define _gkg_qmri_plugin_functors_SingleCompartmentRelaxometryMapper_SingleCompartmentRelaxometryMCMCFunction_h_
#include <gkg-processing-numericalanalysis/MonteCarloMarkovChainFunction.h>
#include <gkg-processing-coordinates/Vector3d.h>
namespace gkg
{
class SingleCompartmentRelaxometryMCMCFunction :
public MonteCarloMarkovChainFunction
{
public:
virtual ~SingleCompartmentRelaxometryMCMCFunction();
double getLikelihoodRatio( const Vector& realMeasurements,
const Vector& currentMeasurements,
const Vector& newMeasurements ) const;
const std::vector< int32_t >& getMeasurementCounts() const;
int32_t getMeasurementCount( int32_t inputVolumeIndex ) const;
const std::vector< double >& getNoiseStandardDeviations() const;
double getNoiseStandardDeviation( int32_t inputVolumeIndex ) const;
const std::vector< double >& getNoiseVariances() const;
double getNoiseVariance( int32_t inputVolumeIndex ) const;
int32_t getParameterCount() const;
void setSite( const Vector3d< int32_t >& site );
virtual void getValuesAt( const Vector& parameters,
Vector& values ) const = 0;
protected:
SingleCompartmentRelaxometryMCMCFunction(
const std::vector< int32_t >& measurementCounts,
const std::vector< double >& noiseStandardDeviations );
std::vector< int32_t > _measurementCounts;
std::vector< double > _noiseStandardDeviations;
int32_t _parameterCount;
std::vector< double > _noiseVariances;
Vector3d< int32_t > _site;
};
}
#endif
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperLoopContext.h>
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMapperGauge.h>
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCFunctionFactory.h>
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryMCMCParameterSet.h>
#include <gkg-qmri-plugin-functors/SingleCompartmentRelaxometryMapper/SingleCompartmentRelaxometryNLPFunction.h>
#include <gkg-processing-container/Volume_i.h>
#include <gkg-processing-numericalanalysis/MonteCarloMarkovChainEstimator.h>
#include <gkg-processing-numericalanalysis/OptimizerConstraint.h>
#include <gkg-processing-algobase/MedianFilter_i.h>
#include <gkg-core-io/Eraser.h>
#include <gkg-core-exception/Exception.h>
#include <iomanip>
gkg::SingleCompartmentRelaxometryMapperLoopContext::
SingleCompartmentRelaxometryMapperLoopContext(
gkg::SingleCompartmentRelaxometryMapperGauge& gauge,
const std::vector< gkg::Volume< float > >& inputVolumes,
const std::vector< gkg::Vector3d< int32_t > >& sites,
const std::string& qMriMappingType,
const std::vector< double >& optimizerParameters,
const std::vector< double >& scalarParameters,
const gkg::RCPointer< gkg::AcquisitionParameterSet >& acquisitionParameterSet,
const int32_t& inputVolumeCount,
const std::vector< int32_t >& measurementCounts,
const int32_t& globalMeasurementCount,
const bool& verbose,
gkg::Volume< float >& protonDensityVolume,
gkg::Volume< float >& relaxationTimeVolume,
const bool& computeFittedMeasurements,
gkg::Volume< float >& fittedMeasurementVolume )
: gkg::LoopContext( &gauge ),
_inputVolumes( inputVolumes ),
_sites( sites ),
_qMriMappingType( qMriMappingType ),
_optimizerParameters( optimizerParameters ),
_scalarParameters( scalarParameters ),
_acquisitionParameterSet( acquisitionParameterSet ),
_inputVolumeCount( inputVolumeCount ),
_measurementCounts( measurementCounts ),
_globalMeasurementCount( globalMeasurementCount ),
_verbose( verbose ),
_protonDensityVolume( protonDensityVolume ),
_relaxationTimeVolume( relaxationTimeVolume ),
_computeFittedMeasurements( computeFittedMeasurements ),
_fittedMeasurementVolume( fittedMeasurementVolume )
{
try
{
////////////////////////////////////////////////////////////////////////////
// getting access to the numerical analysis factory
////////////////////////////////////////////////////////////////////////////
_factory = gkg::NumericalAnalysisSelector::getInstance().
getImplementationFactory();
////////////////////////////////////////////////////////////////////////////
// initializing parameter boundaries, step and fixed information
////////////////////////////////////////////////////////////////////////////
_initialParameters.reallocate( 2 );
_initialParameters( 0 ) = _scalarParameters[ 0 ];
_initialParameters( 1 ) = _scalarParameters[ 1 ];
_lowerParameterBoundaries.reallocate( 2 );
_lowerParameterBoundaries( 0 ) = _scalarParameters[ 2 ];
_lowerParameterBoundaries( 1 ) = _scalarParameters[ 3 ];
_upperParameterBoundaries.reallocate( 2 );
_upperParameterBoundaries( 0 ) = _scalarParameters[ 4 ];
_upperParameterBoundaries( 1 ) = _scalarParameters[ 5 ];
_deltaParameters.reallocate( 2 );
_deltaParameters( 0 ) = _scalarParameters[ 6 ];
_deltaParameters( 1 ) = _scalarParameters[ 7 ];
_isFixedParameters.resize( 2U );
_isFixedParameters[ 0 ] = false;
_isFixedParameters[ 1 ] = false;
////////////////////////////////////////////////////////////////////////////
// preparing NLP constraints and optimizer parameters
////////////////////////////////////////////////////////////////////////////
_constraints.resize( 2U );
_nlpOptimizerParameters.resize( 4U );
_nlpOptimizerParameters[ 0 ] = _optimizerParameters[ 0 ];
_nlpOptimizerParameters[ 1 ] = _optimizerParameters[ 1 ];
int32_t p = 0;
for ( p = 0 ; p < 2; p++ )
{
_constraints[ p ].reset( new gkg::ClosedIntervalOptimizerConstraint(
_lowerParameterBoundaries( p ),
_upperParameterBoundaries( p ) ) );
_nlpOptimizerParameters[ 2 + p ] =
_constraints[ p ]->constrainedDeltaParameter( _deltaParameters( p ) );
}
////////////////////////////////////////////////////////////////////////////
// preparing MCMC optimizer parameters
/////////////////////////////////////////////////////////////////////