/* -*- mona-c++ -*-
* Copyright (c) Leipzig, Madrid 2004 - 2008
* Max-Planck-Institute for Human Cognitive and Brain Science
* Max-Planck-Institute for Evolutionary Anthropology
* BIT, ETSI Telecomunicacion, UPM
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*
*/
// $Id: getmtr.cc 938 2006-07-11 11:57:01Z write1 $
/*! \brief mona-getmtr -- compute magnetization transfer ratio (MTR) images
\par Program description
Computes magnetization transfer ratio (MTR) images, and whole-brain MTR images.
For the measurement a pair of images is used: one without off-resonance saturation pulses,
and one where these pulses are included.
The first image is normally a regular proton-density-weighted image that may be acquired using
a spin-echo or gradient echo sequence. The second image is acquired in an identical way to the
first, except that it includes additional off-resonance RF saturation pulses that affect the
longitudinal magnetisation of protons in the larger molecules, such as cell membranes.
These two images will be referred to as \a M0 and \a Ms respectively.
As it is usually not reasonable to conduct an MTR analysis on the whole image, the analysis comes
in three stages:
- Removel of non-cerebral tissue from the M0 image.
- Calculation of the MTR image.
- Production of the MTR histogram.
The multi-slice \a M0 and \a Ms images need to be given in two separate image files, carried to the
program by the parameters \c --m0-image and \c --ms-image respectively. A region of interest
(e.g. non-cerebral tissue, NAWM, etc) is specified by a binary mask that need to be given as
a seperate image file (\c --mask-image).
For the MTR calculation a decision on an intensity threshold need to be done (\c --threshold).
The MTR calculation wil only be applied to pixels where the intensity is above the threshold.
This prevents the calculation of MTR for pixels that just contain background noise.
The program allows to decide for a regional MTR estimation (<code>--regional-mtr-calc true</code>), or
MTR histogram analysis (<code>--regional-mtr-calc false</code>).
- For regional MTR analysis, the MTR calulation will produce an image with intensities between
0 and 1000. 0 corresponds to 0\% MTR, and 1000 corresponds to 100\% MTR.
This scaling factor is chosen to give a resolution of 0.1\% MTR, which is appropriate for regional
analysis.
- For MTR histogram calculation, it is more approriate to "bin" the histogram into MTR increments of 1%,
and help to smooth the histogram.
Note that the actual histogram calculation need to be done by another program (e.g., \a mona-gethisto), fed by
the output image (\c --out-image ).
\par Magnet Transfer Imaging (General comment)
Magnetisation transfer (MT) imaging investigates macromolecular protons that are essentially invisible with conventional MRI
due to their rapid \f$T_2\f$ relaxation. This 'bound' proton pool can be indirectly imaged by application of a
saturation pre-pulse. Subsequent transfer of saturation from the bound to the free pool occurs via cross-relaxation, spin diffusion,
and chamical exchange. The result of this transfer is to reduce the intensity of the free proton signal.
The size of the MT effect, the MTR, can be quantitatively measured by aquiring two sets of images, with and without
saturation pulse, while all other parameters remain identical.
The MTR of a given region of interest (NAWM, Lesion, etc) is then given by
\f[MTR = ( M_0 - M_s) / M_0\f]
where \f$ M_0 \f$ and \f$M_s\f$ represent the mean pixel intensity of a region without and
with saturation, respectively.
\par Usage
<code>mona-getmtr --usage</code> (or <code> mona-getmtr --help</code>)
\param --m0-image image acquired without saturation pulse
\param --ms-image image acquired with saturation pulse
\param --mask-image a region of interest given as a binary image
\param --threshold threshold value to discriminate background noise
\param --regional-mtr-calc regional MTR estimatatin or pre-processing a MTR histogram
\param --out-image the MTR values staored in an image
\param --verbose some verbose output (min/max values, etc)
\par Example
-# Perform a regional MTR evaluation with two VISTA images
\code
user> mona-getmtr -0 nomtc-image.v -s mtc-image.v -m regional-mask.v -r 1 -o mtr-image.v
\endcode
Note, the image format can be anything covered by the io plugins included in \a libmona.
-# MTR histogram analysis by additionally applying \a mona-gethisto
\code
user> mona-getmtr -0 nomtc-image.v -s mtc-image.v -m regional-mask.v -o - \
moan-histo -i - -o mtr-histogram.dat
\endcode
\par Known bugs
Still none known.
\sa gethisto.cc
\file getmtr.cc
\author M. Tittgemeyer, tittge@cbs.mpg.de, 2004
*/
// $Id: getmtr.cc 938 2006-07-11 11:57:01Z write1 $
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <iostream>
#include <string>
#include <dlfcn.h>
#include <mona.hh>
using namespace std;
using namespace mona;
static const char *revision = "$Revision: 938 $";
#define error_report(level, text) cverb << level << (text) << std::endl;
namespace options {
}
namespace {
} // anonymous namespace
/*! \brief Functor to handle MTR computation
The functor class is initialised by a threshold value to discriminate background noise,
as well as by a scaling factor to the MTR values. The latter is choosen according to the
calculation task, e.g. whether a regional estimation is to perform, of a pre-processing to
further histogram analysis. An optional mask image can be given to assess values in a
specific region of interest.
The class is derived from the image filterbase, that solely defines a result type.
Whenever using in an image filter, the result_type is specified by
<code>typename Filter::result_type</code> that needs to be resolved here.
*/
class CGetMTRatio: public TUnaryImageFilter<C3DImageWrap> {
C3DImageWrap _M_mask;
bool _M_has_mask;
float _M_thresh;
float _M_scaleFactor;
public:
//! \name Constructors
//@{
/*! Default constructor
\param mask a binary mask facilitating a regional analysis
\param threshold threshold value to discriminate background noise
\param regio_calc discriminates the task (regional MTR calc, or histogram pre-pocessing
*/
CGetMTRatio(C3DImageWrap mask, const float& threshold, const bool regio_calc):
_M_mask(mask), _M_has_mask(true), _M_thresh(threshold)
{
_M_scaleFactor = regio_calc ? 1000. : 100.;
}
/*! Alternative constructor
\param threshold threshold value to discriminate background noise
\param regio_calc discriminates the task (regional MTR calc, or histogram pre-pocessing
*/
CGetMTRatio(const float& threshold, const bool regio_calc):
_M_has_mask(false),_M_thresh(threshold)
{
_M_scaleFactor = regio_calc ? 1000. : 100.;
}
//@}
//! \name Operator
//@{
/*! Brace operator
\param dataA any data struture that hold an (STL- style) iterator
\param dataB any data struture that hold an (STL- style) iterator
\returns the computed MTR image (as float float)
*/
template <class Data3D_TYPEA, class Data3D_TYPEB>
typename CGetMTRatio::result_type operator () (Data3D_TYPEA& dataA, Data3D_TYPEB& dataB)
{
// Check whether A and B have the same size
if ( dataA.get_size() != dataB.get_size() ) {
string not_the_same_size = ("Sizes of input data differ ");
error_report(vstream::ml_error, not_the_same_size);
return C3DImageWrap();
};
// Initialize the image supposed to hold the MTR values
C3DFImage *mtr_image = new C3DFImage( dataA.get_size(), dataA.get_attribute_list() );
C3DFImage::iterator __mtr = mtr_image->begin();
// Get data field iterators
typename Data3D_TYPEA::const_iterator __m0 = dataA.begin();
typename Data3D_TYPEB::const_iterator __ms = dataB.begin();
// Get zero value
float Zero = 0.0f;
// some info
cverb << vstream::ml_message << "data scaling by: " << _M_scaleFactor << " with threshold: " << _M_thresh << endl;
// do the job
if ( _M_has_mask ) {
// Check whether A and mask have the same size
if ( dataA.get_size() != _M_mask.get_size() ) {
string not_the_same_size = ("Sizes of input and mask differ ");
error_report(vstream::ml_error, not_the_same_size);
return C3DImageWrap();
};
// get the actual image
C3DBitImage* bit_mask = _M_mask.getC3DBitImage();
// check whether mask is binary
if ( !bit_mask ){
string wrong_type = ("Mask need to be bit type");
error_report(vstream::ml_error, wrong_type);
return C3DImageWrap();
};
typename C3DBitImage::const_iterator __mask = bit_mask->begin();
for ( ; __m0 != dataA.end(); ++__m0, ++__ms, ++__mask, ++__mtr)
*__mtr = ( (*__m0 > _M_thresh) && *__mask ) ?
( (*__m0 - *__ms) * _M_scaleFactor / *__m0 ) : Zero;
} else {
cverb << vstream::ml_message << "NO masking !!" << endl;
for ( ; __m0 != dataA.end(); ++__m0, ++__ms, ++__mtr)
*__mtr = ( *__m0 > _M_thresh ) ?
( (*__m0 - *__ms) * _M_scaleFactor / *__m0 ) : Zero;
};
// get result
return C3DImageWrap( mtr_image );
}
//@}
}; // CGetMTRatio
int main( int argc, const char *argv[] )
{
string m0_filename;
string ms_filename;
string mask_filename;
string mtr_filename = "-";
float threshold = 100.0f;
bool regio_calc = false;
try {
popt::COptions opts;
opts.push_back(popt::option( m0_filename, "m0-image", 'i',
"image acquired without saturation pulse", NULL ));
opts.push_back(popt::option( ms_filename, "ms-image", 's',
"image acquired with additional saturation pulse", NULL ));
opts.push_back(popt::option( mtr_filename, "mtr-image", 'o',
"image holding the MTR values", "-" ));
opts.push_back(popt::option( mask_filename, "mask-image", 'm',
"a region of interest given as binary image mask", NULL ));
opts.push_back(popt::option( threshold, "threshold", 't',
"threshold value to dicriminate background noise", "100" ));
opts.push_back(popt::option( regio_calc, "regional-mtr-calc", 'r',
"to perform a regional MTR calcultion (false, MTR-histogram pre-processing" ));
vector<string> non_options;
popt::parse_options(argc, argv, opts, non_options);
if ( non_options.size() > 0 )
throw invalid_argument("unknown options");
if ( m0_filename.empty() )
throw mona_runtime_error("'--m0-image' ('i') option required");
if ( ms_filename.empty() )
throw mona_runtime_error("'--ms-image' ('m') option required");
C3DImageIOPluginHandler imageio;
// initialize the functor for MTR- calculation
auto_ptr<CGetMTRatio> mtrRatio;
if ( mask_filename.empty() ) {
mtrRatio.reset( new CGetMTRatio(threshold, regio_calc) );
} else {
// read mask
auto_ptr<C3DImageList> mask_image_list(imageio.load(mask_filename));
if (!mask_image_list.get() || (mask_image_list->size() == 0)) {
string not_found = ("No supported data found in ") + mask_filename;
throw mona_runtime_error(not_found);
};
mtrRatio.reset( new CGetMTRatio(*mask_image_list->begin(), threshold, regio_calc) );
}
// read M0
auto_ptr<C3DImageList> m0_image_list(imageio.load(m0_filename));
if (!m0_image_list.get() || !m0_image_list->size()) {
string not_found = ("No supported data found in ") + m0_filename;
throw mona_runtime_error(not_found);
}
// read Ms
auto_ptr<C3DImageList> ms_image_list(imageio.load(ms_filename));
if (!ms_image_list.get() || !ms_image_list->size()) {
string not_found = ("No supported data found in ") + ms_filename;
throw mona_runtime_error(not_found);
}
// get the MTR
C3DImageWrap mtr_image( image_filter( *mtrRatio, *m0_image_list->begin(), *ms_image_list->begin() ));
// save to file
C3DImageList out_list;
out_list.push_back( mtr_image );
// write history
CHistory::instance().append(argv[0], revision, opts);
if ( !imageio.save(m0_image_list->get_sourceformat(), out_list, mtr_filename) ){
string not_save = ("unable to save result to ") +
mtr_filename;
throw mona_runtime_error(not_save);
};
return EXIT_SUCCESS;
}
catch (const mona_runtime_error& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const mona_fatal_error& e){
cerr << argv[0] << " fatal: " << e.what() << endl;
}
catch (const mona_exception& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const invalid_argument &e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (const exception& e){
cerr << argv[0] << " error: " << e.what() << endl;
}
catch (...){
cerr << argv[0] << " unknown exception" << endl;
}
return EXIT_FAILURE;
}
/*
$Log$
Revision 1.20 2005/06/29 13:32:02 wollny
move to libmona-version 0.7
Revision 1.2 2005/06/02 13:33:23 gerddie
adapt code to new plugin handling
Revision 1.1.1.1 2005/03/17 13:48:32 gerddie
initial checkin
Revision 1.19 2005/01/27 17:46:51 tittge
remove bug that caused segfault when not using the masking option
Revision 1.18 2005/01/19 11:10:07 wollny
adapted to new verbose handling of libmona
Revision 1.17 2005/01/18 12:54:22 wollny
Changes in option handling made it neccessary to update all programs. The
default vale of an option is now given as the last !string! vale of that
option.
Revision 1.16 2004/11/22 10:23:06 tittge
adapt functors
Revision 1.15 2004/11/01 16:24:54 wollny
move from old error-report to cverb
Revision 1.14 2004/10/22 15:03:04 tittge
remove some memory problems (by consitently using auto_ptr's)
Revision 1.13 2004/10/15 14:24:45 tittge
added a crisp (maximum likelyhood) segmentation
Revision 1.12 2004/08/23 12:00:00 wollny
adapt to warp-image
Revision 1.11 2004/08/11 07:26:45 wollny
added support for multiple images
Revision 1.10 2004/06/03 09:57:45 wollny
Changed (hopefully) all instancable class names to Cxxxxx
Revision 1.9 2004/05/13 14:29:26 wollny
adapt to new message report
Revision 1.8 2004/04/05 08:44:04 tittge
make dot tool compatible
Revision 1.7 2004/03/24 08:19:00 tittge
change exit value (for Makefile usage)
Revision 1.6 2004/03/23 16:03:17 tittge
some bugs removed
Revision 1.5 2004/03/23 13:16:47 tittge
functors re-defined
Revision 1.4 2004/03/22 08:44:28 tittge
extensive re-programmimng on histogram
Revision 1.3 2004/03/19 14:08:44 tittge
changes in coding style
Revision 1.2 2004/03/17 14:35:01 tittge
program redesiged
Revision 1.1 2004/03/15 08:04:53 tittge
program for calculating MTR
Revision 1.1 2004/03/12 09:31:54 tittge
modified program structure
Revision 1.2 2004/02/23 15:45:04 tittge
update to libmona-0.1 functionality
*/