diff --git a/UKFTractography/UKFTractography.cxx b/UKFTractography/UKFTractography.cxx index 4988406..7a5f701 100644 --- a/UKFTractography/UKFTractography.cxx +++ b/UKFTractography/UKFTractography.cxx @@ -63,9 +63,14 @@ int ModuleEntryPoint(int argc, char **argv) try { if (tract->LoadFiles(ukf_settings.dwiFile, - ukf_settings.seedsFile, ukf_settings.maskFile, - normalizedDWIData, outputNormalizedDWIData) == EXIT_FAILURE) + normalizedDWIData, + outputNormalizedDWIData == EXIT_FAILURE, + ukf_settings.seedsFile, + ukf_settings.stopFile, + ukf_settings.wmFile, + ukf_settings.gmFile, + ukf_settings.csfFile)) { itkGenericExceptionMacro(<< "::LoadFiles failed with unknown error."); } diff --git a/UKFTractography/UKFTractography.xml b/UKFTractography/UKFTractography.xml index 0acbad4..1c08fa4 100644 --- a/UKFTractography/UKFTractography.xml +++ b/UKFTractography/UKFTractography.xml @@ -14,7 +14,6 @@ Yogesh Rathi, Stefan Lienhard, Yinpeng Li, Martin Styner, Ipek Oguz, Yundi Shi, Christian Baumgartner, Kent Williams, Hans Johnson, Peter Savadjiev, Carl-Fredrik Westin, Lauren O'Donnell, Jessica Lee. - Input/output parameters @@ -27,26 +26,10 @@ Input diffusion weighted (DWI) volume - - seedsFile - seedsFile - - input - Seeds for diffusion. If not specified, full brain tractography will be performed, and the algorithm will start from every voxel in the brain mask where the Generalized Anisotropy is bigger than 0.18 - - - - labels - labels - - A vector of the ROI labels to be used. There are the voxel values where tractography should be seeded. - 1 - - maskFile maskFile - + input Brain mask for diffusion tractography. Tracking will only be performed inside this mask. @@ -62,30 +45,17 @@ - - Basic Parameters - - - seedsPerVoxel - seedsPerVoxel - - Tractography parameter used in all models. Each seed generates a fiber, thus using more seeds generates more fibers. In general use 1 or 2 seeds, and for a more thorough result use 5 or 10 (depending on your machine this may take up to 2 days to run). Default: 1. Range: 0-50. - - 1 - - 0 - 50 - 1 - - + + Options for seeding tractography. Only one of the three provided options will be used. Option 1 (Minimum Seed FA) is the default. + seedingThreshold seedingThreshold - - Tractography parameter used in all models. Seed points whose fractional anisotropy (FA) or mean signal are below this value are excluded. Default: 0.18. Range: 0-1. + + (Seeding Option 1) Tractography parameter used in all models. Seed points whose fractional anisotropy (FA) are below this value are excluded. This seeding option is default and will not be used if Seeding Option 2 or 3 is used. Default: 0.18. Range: 0-1. 0.18 0 @@ -94,13 +64,57 @@ + + wmFile + wmFile + + input + (Seeding Option 2) A probabilistic segmentation map of the White Matter (WM). The values in the map should be between 0 and 1. Voxels with a probability over WM Probability Threshold (--wmProbThreshold) will be used for seeding. This seeding option will not be used if Seeding Option 3 is used. + + + + wmProbThreshold + wmProbThreshold + + (Seeding Option 2) When a WM Segmentation Map (--wmFile) is provided, tracking will be seeded in voxels with values over this threshold. Default: 0.3. Range: 0-1. + 0.3 + + 0 + 1 + 0.01 + + + + + seedsFile + seedsFile + + input + (Seeding Option 3) Voxels in this map with a label defined in ROI Labels (--seedLabels) for Seeding will be used for seeding. If this option is used, Seeding Options 1 and 2 will not be used. + + + + seedLabels + seedLabels + + (Seeding Option 3) A list of the ROI labels to be used when Seeding Label Map (--seedsFile) is provided. There are the voxel values where tractography should be seeded. + 1 + + + + + + + + Options for stopping tractography. Only one of the three provided options will be used. Option 1 (Terminating Mean Signal) is the default. + stoppingFA stoppingFA - - Tractography parameter used in tensor model. Tractography will stop when the fractional anisotropy (FA) of the tensor being tracked is less than this value. Note: make sure to also decrease the GA to track through lower anisotropy areas. This parameter is used only in tensor models. Default: 0.15. Range: 0-1. + + (Stopping Option 1) Tractography parameter used only in tensor models. Tractography will stop when the fractional anisotropy (FA) of the tensor being tracked is less than this value. Note: make sure to also decrease the Terminating Mean Signal (--stoppingThreshold) to track through lower anisotropy areas. This option will not be used if Stopping Option 2 or 3 is provided. Default: 0.15. Range: 0-1. 0.15 0 @@ -112,8 +126,8 @@ stoppingThreshold stoppingThreshold - - Tractography parameter used in all models. Tractography will stop when the mean signal is below this value. Default: 0.1. Range: 0-1. + + (Stopping Option 1) Tractography parameter used by default in all models. Tractography will stop when the mean signal is below this value. This option will not be used if Stopping Option 2 or 3 is provided. Default: 0.1. Range: 0-1. 0.1 0 @@ -122,6 +136,84 @@ + + gmFile + gmFile + + input + (Stopping Option 2) A probabilistic segmentation map of the Gray Matter (GM). The values in the map should be between 0 and 1. Tracking will stop in the voxels with a probability over GM Probability Threshold (--gmProbThreshold). This option will not be used if Stopping Option 3 is used. + + + + gmProbThreshold + gmProbThreshold + + (Stopping Option 2) When a GM Segmentation Map (--gmFile) is provided, tracking will stop in voxels with values over this threshold. Default: 0.99. Range: 0-1. + 0.99 + + 0 + 1 + 0.01 + + + + + csfFile + csfFile + + input + (Stopping Option 2) A probabilistic segmentation map of the Cerebrospinal Fluid (CSF). The values in the map should be between 0 and 1. Tracking will stop in voxels with a probability over CSF Probability Threshold (--csfProbThreshold). This option will not be used if Stopping Option 3 is used. + + + + csfProbThreshold + csfProbThreshold + + (Stopping Option 2) When a CSF Segmentation Map (--csfFile) is provided, tracking will stop in voxels with values over this threshold. Default: 0.5. Range: 0-1. + 0.5 + + 0 + 1 + 0.01 + + + + + stopFile + stopFile + + input + (Stopping Option 3) Label map that defines where tracking should stop. Voxels in this map with a label listed in ROI Labels for Stopping (--stopLabels) will be used for stopping. If this option is provided, Stopping Options 1 and 2 will not be used. + + + + stopLabels + stopLabels + + (Stopping Option 3) A list of the ROI labels to be used when Stopping Label Map (--stopFile) is provided. There are the voxel values where tractography should stop. + 1 + + + + + + + Basic Parameters + + + seedsPerVoxel + seedsPerVoxel + + Tractography parameter used in all models. Each seed generates a fiber, thus using more seeds generates more fibers. In general use 1 or 2 seeds, and for a more thorough result use 5 or 10 (depending on your machine this may take up to 2 days to run). Default: 1. Range: 0-50. + + 1 + + 0 + 50 + 1 + + + numThreads numThreads @@ -211,7 +303,7 @@ freeWater freeWater - Adds a term for free water diffusion to the model. The free water model is a tensor with all 3 eigenvalues equal to the diffusivity of free water (0.003). To output the free water fraction, make sure to use the "save free water" parameter. + Adds a term for free water diffusion to the model. The free water model is a tensor with all 3 eigenvalues equal to the diffusivity of free water (0.003). To output the free water fraction, make sure to use the "save free water" parameter. false diff --git a/ukf/ISignalData.h b/ukf/ISignalData.h index ffeabf3..8e68bd8 100644 --- a/ukf/ISignalData.h +++ b/ukf/ISignalData.h @@ -50,8 +50,14 @@ class ISignalData /** Checks if a certian position is still within the brain mask. */ virtual ukfPrecisionType ScalarMaskValue(const vec3_t& pos) const = 0; + /** Gets stopping mask value at a certain position, from stop label map, GM segmentation, and CSF segmentation */ + virtual ukfPrecisionType ScalarStopValue(const std::vector& labels, const ukfPrecisionType gm_prob_threshold, const ukfPrecisionType csf_prob_threshold, const vec3_t& pos) const = 0; + + /** Check if GM or CSF masks are provided for stopping */ + virtual ukfPrecisionType isGMCSFProvided() const = 0; + /** Get all the seed points. */ - virtual void GetSeeds(const std::vector& labels, stdVec_t& seeds) const = 0; + virtual void GetSeeds(const std::vector& labels, const ukfPrecisionType prob_threshold, stdVec_t& seeds) const = 0; /** Returns the gradients. */ virtual const stdVec_t & gradients() const = 0; @@ -72,8 +78,10 @@ class ISignalData * * Loads all the data necessary to perform tractography */ - virtual bool LoadData(const std::string& data_file, const std::string& seed_file, const std::string& mask_file, - const bool normalizedDWIData, const bool outputNormalizedDWIData) = 0; + virtual bool LoadData(const std::string& data_file, const std::string& mask_file, + const bool normalizedDWIData, const bool outputNormalizedDWIData, + const std::string& seed_file, const std::string& stop_file, + const std::string& wm_file, const std::string& gm_file, const std::string& csf_file) = 0; /** Returns the dimensions of the image */ virtual vec3_t dim() const = 0; diff --git a/ukf/NrrdData.cc b/ukf/NrrdData.cc index fb623a4..842d633 100644 --- a/ukf/NrrdData.cc +++ b/ukf/NrrdData.cc @@ -12,7 +12,7 @@ NrrdData::NrrdData(ukfPrecisionType sigma_signal, ukfPrecisionType sigma_mask) : ISignalData(sigma_signal, sigma_mask), - _data(NULL), _seed_data(NULL), _mask_data(NULL), _data_nrrd(NULL) + _data(NULL), _seed_data(NULL), _stop_data(NULL), _mask_data(NULL), _wm_data(NULL), _gm_data(NULL), _csf_data(NULL), _data_nrrd(NULL) { } @@ -26,6 +26,22 @@ NrrdData::~NrrdData() { nrrdNuke(_seed_nrrd); } + if( _stop_data ) + { + nrrdNuke(_stop_nrrd); + } + if( _wm_data ) + { + nrrdNuke(_wm_nrrd); + } + if( _gm_data ) + { + nrrdNuke(_gm_nrrd); + } + if( _csf_data ) + { + nrrdNuke(_csf_nrrd); + } if( _mask_data ) { nrrdNuke(_mask_nrrd); @@ -125,7 +141,7 @@ ukfPrecisionType NrrdData::ScalarMaskValue(const vec3_t& pos) const const int nz = static_cast(_dim[2]); unsigned int index; - ukfPrecisionType value; + ukfPrecisionType value; const int x = static_cast(round(pos[0])); const int y = static_cast(round(pos[1])); @@ -154,13 +170,187 @@ ukfPrecisionType NrrdData::ScalarMaskValue(const vec3_t& pos) const } break; default: - std::cout << "Unsupported data type for seed file!" << std::endl; + std::cout << "Unsupported data type for mask file!" << std::endl; throw; } return value; } +ukfPrecisionType NrrdData::isGMCSFProvided() const +{ + if ( _gm_data || _csf_data || _stop_data) + { + return ukfOne; + } + else + { + return ukfZero; + } +} + +ukfPrecisionType NrrdData::ScalarStopValue(const std::vector& stopLabels, const ukfPrecisionType gm_prob_threshold, const ukfPrecisionType csf_prob_threshold, const vec3_t& pos) const +{ + const int nx = static_cast(_dim[0]); + const int ny = static_cast(_dim[1]); + const int nz = static_cast(_dim[2]); + + unsigned int index; + ukfPrecisionType is_stopping_voxel = ukfZero; + + const int x = static_cast(round(pos[0])); + const int y = static_cast(round(pos[1])); + const int z = static_cast(round(pos[2])); + + if( (x < 0 || nx <= x) || + (y < 0 || ny <= y) || + (z < 0 || nz <= z) ) + { + return is_stopping_voxel; + } + + index = nz * ny * x + nz * y + z; + + if ( _stop_data ) + { + + size_t nx_stop = _stop_nrrd->axis[2].size; + size_t ny_stop = _stop_nrrd->axis[1].size; + size_t nz_stop = _stop_nrrd->axis[0].size; + assert(_stop_data); + + if ( !(nx_stop == _dim[0] && ny_stop == _dim[1] && nz_stop == _dim[2]) ) + { + itkGenericExceptionMacro(<< "Stopping Labelmap ROI volume dimensions DO NOT match DWI dimensions"); + } + + // signed or unsigned doesn't make a difference since masks don't contain any negative values + int stop_value = 0; + switch( _stop_data_type ) + { + case 0: + { + stop_value = 0; + } + break; + case 2: + { + stop_value = static_cast(_stop_data)[index]; + } + break; + case 3: + { + stop_value = static_cast(_stop_data)[index]; + } + break; + case 5: + { + stop_value = static_cast(_stop_data)[index]; + } + break; + default: + std::cout << "Unsupported data type for stop file!" << std::endl; + throw; + } + + std::vector::const_iterator cit; + for( cit = stopLabels.begin(); cit != stopLabels.end(); ++cit ) + { + if( *cit == stop_value ) + { + is_stopping_voxel = ukfOne; + } + } + + } + else // If stopping label map is provided, GM and CSF will not be used. + { + if ( _gm_data ) + { + size_t nx_gm = _gm_nrrd->axis[2].size; + size_t ny_gm = _gm_nrrd->axis[1].size; + size_t nz_gm = _gm_nrrd->axis[0].size; + assert(_gm_data); + + if ( !(nx_gm == _dim[0] && ny_gm == _dim[1] && nz_gm == _dim[2]) ) + { + itkGenericExceptionMacro(<< "GM segmentation volume dimensions DO NOT match DWI dimensions"); + } + + float gm_prob = 0.0; + switch( _gm_data_type ) + { + case 0: + { + gm_prob = 0.0; + } + break; + case 9: + { + gm_prob = static_cast(_gm_data)[index]; + } + break; + case 10: + { + gm_prob = static_cast(_gm_data)[index]; + } + break; + default: + std::cout << "Unsupported data type for GM segmentation file!" << std::endl; + throw; + } + + if (gm_prob > gm_prob_threshold) + { + is_stopping_voxel = ukfOne; + } + } + + if ( _csf_data ) + { + size_t nx_csf = _csf_nrrd->axis[2].size; + size_t ny_csf = _csf_nrrd->axis[1].size; + size_t nz_csf = _csf_nrrd->axis[0].size; + assert(_csf_data); + + if ( !(nx_csf == _dim[0] && ny_csf == _dim[1] && nz_csf == _dim[2]) ) + { + itkGenericExceptionMacro(<< "CSF segmentation volume dimensions DO NOT match DWI dimensions"); + } + + float csf_prob = 0.0; + switch( _csf_data_type ) + { + case 0: + { + csf_prob = 0.0; + } + break; + case 9: + { + csf_prob = static_cast(_csf_data)[index]; + } + break; + case 10: + { + csf_prob = static_cast(_csf_data)[index]; + } + break; + default: + std::cout << "Unsupported data type for GM segmentation file!" << std::endl; + throw; + } + + if (csf_prob > csf_prob_threshold) + { + is_stopping_voxel = ukfOne; + } + } + } + + return is_stopping_voxel; +} + ukfPrecisionType NrrdData::Interp3ScalarMask(const vec3_t& pos) const { const int nx = static_cast(_dim[0]); @@ -235,74 +425,133 @@ ukfPrecisionType NrrdData::Interp3ScalarMask(const vec3_t& pos) const return s / w_sum; } -void NrrdData::GetSeeds(const std::vector& labels, +void NrrdData::GetSeeds(const std::vector& labels, const ukfPrecisionType wm_prob_threshold, stdVec_t& seeds) const { + + const int nx = static_cast(_dim[0]); + const int ny = static_cast(_dim[1]); + const int nz = static_cast(_dim[2]); + if( _seed_data ) { - assert(seeds.size() == 0); - std::vector::const_iterator cit; + std::cout << "Seeding using Option 3, from a seeding label map." << std::endl; + size_t nx_seed = _seed_nrrd->axis[2].size; + size_t ny_seed = _seed_nrrd->axis[1].size; + size_t nz_seed = _seed_nrrd->axis[0].size; + assert(_seed_data); + + if ( !(nx_seed == _dim[0] && ny_seed == _dim[1] && nz_seed == _dim[2]) ) + { + itkGenericExceptionMacro(<< "Seeding Labelmap ROI volume dimensions DO NOT match DWI dimensions"); + } + } + else if( _wm_data ) + { + std::cout << "Seeding using Option 2, from a WM segmentation map." << std::endl;; + size_t nx_wm = _wm_nrrd->axis[2].size; + size_t ny_wm = _wm_nrrd->axis[1].size; + size_t nz_wm = _wm_nrrd->axis[0].size; + assert(_wm_data); - // Go through the volume. - size_t nx = _seed_nrrd->axis[2].size; - size_t ny = _seed_nrrd->axis[1].size; - size_t nz = _seed_nrrd->axis[0].size; - assert(_seed_data); + if ( !(nx_wm == _dim[0] && ny_wm == _dim[1] && nz_wm == _dim[2]) ) + { + itkGenericExceptionMacro(<< "WM segmentation volume dimensions DO NOT match DWI dimensions"); + } + } - if ( !(nx == _dim[0] && ny == _dim[1] && nz == _dim[2]) ) - { - itkGenericExceptionMacro(<< "Labelmap ROI volume dimensions DO NOT match DWI dimensions"); - } + assert(seeds.size() == 0); // Make sure seeds is empty - for( size_t i = 0; i < nx; ++i ) + for( int i = 0; i < nx; ++i ) + { + for( int j = 0; j < ny; ++j ) { - for( size_t j = 0; j < ny; ++j ) + for( int k = 0; k < nz; ++k ) { - for( size_t k = 0; k < nz; ++k ) - { - for( cit = labels.begin(); cit != labels.end(); ++cit ) + int is_seeding_voxel = 0; + size_t index = ny * nz * i + nz * j + k; + + if( _seed_data ) { - int value = 0; - size_t index = ny * nz * i + nz * j + k; + std::vector::const_iterator cit; + for( cit = labels.begin(); cit != labels.end(); ++cit ) + { + int value = 0; - switch( _seed_data_type ) + switch( _seed_data_type ) + { + case 2: + { + value = static_cast(_seed_data)[index]; + } + break; + case 3: + { + value = static_cast(_seed_data)[index]; + } + break; + case 5: + { + value = static_cast(_seed_data)[index]; + } + break; + default: + std::cout << "Unsupported data type for seed file!" << std::endl; + assert(false); + } + if( *cit == value ) + { + is_seeding_voxel = 1; + } + } + } + else if ( _wm_data ) + { + + float wm_prob = 0.0; + + switch( _wm_data_type ) { - case 2: + case 0: { - value = static_cast(_seed_data)[index]; + wm_prob = 0.0; } break; - case 3: + case 9: { - value = static_cast(_seed_data)[index]; + wm_prob = static_cast(_wm_data)[index]; } break; - case 5: + case 10: { - value = static_cast(_seed_data)[index]; + wm_prob = static_cast(_wm_data)[index]; } break; default: - std::cout << "Unsupported data type for seed file!" << std::endl; - assert(false); + std::cout << "Unsupported data type for WM segmentation file!" << std::endl; + throw; } - if( *cit == value ) + + if (wm_prob > wm_prob_threshold) // TODO: may support a user provided threshold. { - seeds.push_back(vec3_t(i, j, k) ); + is_seeding_voxel = 1; } } - } + + if( is_seeding_voxel == 1 ) + { + seeds.push_back(vec3_t(i, j, k) ); + } } } } - else - { - std::cout << "No seed data available." << std::endl; - } -} + + } -bool NrrdData::SetData(Nrrd* data_nrrd, Nrrd* mask_nrrd, Nrrd* seed_nrrd, - bool normalizedDWIData) +bool NrrdData::SetData(Nrrd* data_nrrd, Nrrd* mask_nrrd, + bool normalizedDWIData, + Nrrd* seed_nrrd, Nrrd* stop_nrrd, + Nrrd* wm_nrrd, Nrrd* gm_nrrd, Nrrd* csf_nrrd) { //_data_nrrd = (Nrrd*)data; //_seed_data = (Nrrd*)seed; @@ -313,14 +562,6 @@ bool NrrdData::SetData(Nrrd* data_nrrd, Nrrd* mask_nrrd, Nrrd* seed_nrrd, return true; } - if (seed_nrrd) - { - this->_seed_nrrd = seed_nrrd; - this->_seed_data = seed_nrrd->data; - this->_seed_data_type = seed_nrrd->type; - assert(_seed_data_type == 2 || _seed_data_type == 3 || _seed_data_type == 5); - } - if( mask_nrrd->type == 1 || mask_nrrd->type == 2 ) { this->_mask_num_bytes = 1; @@ -340,17 +581,60 @@ bool NrrdData::SetData(Nrrd* data_nrrd, Nrrd* mask_nrrd, Nrrd* seed_nrrd, this->_mask_nrrd = mask_nrrd; this->_mask_data = mask_nrrd->data; + if (seed_nrrd) + { + this->_seed_nrrd = seed_nrrd; + this->_seed_data = seed_nrrd->data; + this->_seed_data_type = seed_nrrd->type; + assert(_seed_data_type == 2 || _seed_data_type == 3 || _seed_data_type == 5); + } + + if (stop_nrrd) + { + this->_stop_nrrd = stop_nrrd; + this->_stop_data = stop_nrrd->data; + this->_stop_data_type = stop_nrrd->type; + assert(_stop_data_type == 2 || _stop_data_type == 3 || _stop_data_type == 5); + } + + if (wm_nrrd) + { + this->_wm_nrrd = wm_nrrd; + this->_wm_data = wm_nrrd->data; + this->_wm_data_type = wm_nrrd->type; + assert(_wm_data_type == 9 || _wm_data_type == 10); + } + + if (gm_nrrd) + { + this->_gm_nrrd = gm_nrrd; + this->_gm_data = gm_nrrd->data; + this->_gm_data_type = gm_nrrd->type; + assert(_gm_data_type == 9 || _gm_data_type == 10); + } + + if (csf_nrrd) + { + this->_csf_nrrd = csf_nrrd; + this->_csf_data = csf_nrrd->data; + this->_csf_data_type = csf_nrrd->type; + assert(_csf_data_type == 9 || _csf_data_type == 10); + } + return false; } bool NrrdData::LoadData(const std::string& data_file, - const std::string& seed_file, const std::string& mask_file, const bool normalizedDWIData, - const bool outputNormalizedDWIData - ) + const bool outputNormalizedDWIData, + const std::string& seed_file, + const std::string& stop_file, + const std::string& wm_file, + const std::string& gm_file, + const std::string& csf_file) { - if( _data || _seed_data || _mask_data ) + if( _data || _seed_data || _stop_data || _mask_data || _wm_data || _gm_data || _csf_data ) { std::cout << "There is already some data!" << std::endl; return true; @@ -365,6 +649,16 @@ bool NrrdData::LoadData(const std::string& data_file, return true; } + Nrrd* mask_nrrd = nrrdNew(); + // Load mask + if( nrrdLoad(mask_nrrd, mask_file.c_str(), NULL) ) + { + char *err = biffGetDone(NRRD); + std::cout << "Trouble reading " << mask_file << ": " << err << std::endl; + free( err ); + return true; + } + if( outputNormalizedDWIData ) { std::string normalizedDataFile = data_file.substr(0, data_file.find_last_of('.') ); @@ -380,9 +674,8 @@ bool NrrdData::LoadData(const std::string& data_file, } } - Nrrd* seed_nrrd = nrrdNew(); - // Load seeds + // Load seeding map if( !seed_file.empty() ) { if( nrrdLoad(seed_nrrd, seed_file.c_str(), NULL) ) @@ -394,17 +687,59 @@ bool NrrdData::LoadData(const std::string& data_file, } } - Nrrd* mask_nrrd = nrrdNew(); - // Load mask - if( nrrdLoad(mask_nrrd, mask_file.c_str(), NULL) ) + Nrrd* stop_nrrd = nrrdNew(); + // Load stopping map + if( !stop_file.empty() ) { - char *err = biffGetDone(NRRD); - std::cout << "Trouble reading " << mask_file << ": " << err << std::endl; - free( err ); - return true; + if( nrrdLoad(stop_nrrd, stop_file.c_str(), NULL) ) + { + char *err = biffGetDone(NRRD); + std::cout << "Trouble reading " << stop_file << ": " << err << std::endl; + free( err ); + return true; + } + } + + Nrrd* wm_nrrd = nrrdNew(); + // Load WM segmentation + if( !wm_file.empty() ) + { + if( nrrdLoad(wm_nrrd, wm_file.c_str(), NULL) ) + { + char *err = biffGetDone(NRRD); + std::cout << "Trouble reading " << wm_file << ": " << err << std::endl; + free( err ); + return true; + } + } + + Nrrd* gm_nrrd = nrrdNew(); + // Load GM segmentation + if( !gm_file.empty() ) + { + if( nrrdLoad(gm_nrrd, gm_file.c_str(), NULL) ) + { + char *err = biffGetDone(NRRD); + std::cout << "Trouble reading " << gm_file << ": " << err << std::endl; + free( err ); + return true; + } + } + + Nrrd* csf_nrrd = nrrdNew(); + // Load CSF segmentation + if( !csf_file.empty() ) + { + if( nrrdLoad(csf_nrrd, csf_file.c_str(), NULL) ) + { + char *err = biffGetDone(NRRD); + std::cout << "Trouble reading " << csf_file << ": " << err << std::endl; + free( err ); + return true; + } } - bool status = SetData(input_nrrd, mask_nrrd, seed_nrrd, normalizedDWIData); + bool status = SetData(input_nrrd, mask_nrrd, normalizedDWIData, seed_nrrd, stop_nrrd, wm_nrrd, gm_nrrd, csf_nrrd); return status; } @@ -438,7 +773,7 @@ bool NrrdData::LoadSignal(Nrrd* input_nrrd, const bool normalizedDWIData) } assert(size == 2); - ukfPrecisionType bValue = ukfZero; + ukfPrecisionType bValue = ukfZero; assert(_gradients.size() == 0); // Read key value pairs. diff --git a/ukf/NrrdData.h b/ukf/NrrdData.h index 2693403..7406e7c 100644 --- a/ukf/NrrdData.h +++ b/ukf/NrrdData.h @@ -39,15 +39,21 @@ class NrrdData : public ISignalData /** Gets brain mask value at a certain position */ virtual ukfPrecisionType ScalarMaskValue(const vec3_t& pos) const; + /** Gets stopping mask value at a certain position, from stop label map, GM segmentation, and CSF segmentation */ + virtual ukfPrecisionType ScalarStopValue(const std::vector& labels, const ukfPrecisionType gm_prob_threshold, const ukfPrecisionType csf_prob_threshold, const vec3_t& pos) const; + + /** Check if GM or CSF masks are provided for stopping */ + virtual ukfPrecisionType isGMCSFProvided() const; + /** - * \brief Get the seed points from the nrrd file + * \brief Get the seed points from the nrrd files: seeding label map and WM segmentation * * Takes care of different seed data types by type casting * * \param[in] labels a vector of labels that define the seed region * \param[out] seeds a vector containing the positions in ijk-space of the seeds */ - virtual void GetSeeds(const std::vector& labels, stdVec_t& seeds) const; + virtual void GetSeeds(const std::vector& labels, const ukfPrecisionType wm_prob_threshold, stdVec_t& seeds) const; /** returns the gradients of the diffusion image */ virtual const stdVec_t & gradients() const @@ -84,11 +90,13 @@ class NrrdData : public ISignalData * * Loads all the data necessary to perform tractography */ - virtual bool LoadData(const std::string& data_file, const std::string& seed_file, const std::string& mask_file, - const bool normalizedDWIData, const bool outputNormalizedDWIData); + virtual bool LoadData(const std::string& data_file, const std::string& mask_file, + const bool normalizedDWIData, const bool outputNormalizedDWIData, + const std::string& seed_file, const std::string& stop_file, + const std::string& wm_file, const std::string& gm_file, const std::string& csf_file); - virtual bool SetData(Nrrd* data, Nrrd* seed, Nrrd* mask, bool normalizedDWIData); + virtual bool SetData(Nrrd* data, Nrrd* mask, bool normalizedDWIData, Nrrd* seed, Nrrd* stop, Nrrd* wm, Nrrd* gm, Nrrd* csf); /** Returns the dimensions of the signal in each directions as a vector */ vec3_t dim() const @@ -119,10 +127,26 @@ class NrrdData : public ISignalData void *_seed_data; /** seed type is needed for correct casting type */ int _seed_data_type; + /** pointer to stop data, is casted at runtime */ + void *_stop_data; + /** stop type is needed for correct casting type */ + int _stop_data_type; /** pointer to mask data, is casted at runtime */ void *_mask_data; /** number of bytes of the mask is needed for casting */ int _mask_num_bytes; + /** pointer to WM data, is casted at runtime */ + void *_wm_data; + /** WM type is needed for correct casting type */ + int _wm_data_type; + /** pointer to GM data, is casted at runtime */ + void *_gm_data; + /** GM type is needed for correct casting type */ + int _gm_data_type; + /** pointer to CSF data, is casted at runtime */ + void *_csf_data; + /** CSF type is needed for correct casting type */ + int _csf_data_type; /** the actual diffusion data in Nrrd type */ Nrrd *_data_nrrd; @@ -130,6 +154,14 @@ class NrrdData : public ISignalData Nrrd *_seed_nrrd; /** The actual mask data */ Nrrd *_mask_nrrd; + /** The actual stop data */ + Nrrd *_stop_nrrd; + /** The actual WM data */ + Nrrd *_wm_nrrd; + /** The actual GM data */ + Nrrd *_gm_nrrd; + /** The actual CSF data */ + Nrrd *_csf_nrrd; }; #endif // NRRDDATA_H_ diff --git a/ukf/cli.cc b/ukf/cli.cc index cc807e6..bf2b0cc 100644 --- a/ukf/cli.cc +++ b/ukf/cli.cc @@ -125,8 +125,8 @@ int ukf_parse_cli(int argc, char** argv, UKFSettings& s) l_maxBranchingAngle = 0.0; } - if (labels.size() == 0) { - labels.push_back(1) ; //Default to use label 1 + if (seedLabels.size() == 0) { + seedLabels.push_back(1) ; //Default to use label 1 } if (l_stoppingFA == 0.15) { @@ -274,6 +274,9 @@ int ukf_parse_cli(int argc, char** argv, UKFSettings& s) s.transform_position = true; // TODO hard-coded :/ s.store_glyphs = storeGlyphs; s.branches_only = false; // TODO hard-coded :/ + s.wm_prob_threshold = wmProbThreshold; + s.gm_prob_threshold = gmProbThreshold; + s.csf_prob_threshold = csfProbThreshold; s.fa_min = l_stoppingFA; s.mean_signal_min = l_stoppingThreshold; s.seeding_threshold = l_seedingThreshold; @@ -287,7 +290,8 @@ int ukf_parse_cli(int argc, char** argv, UKFSettings& s) s.stepLength = l_stepLength; s.recordLength = l_recordLength; s.maxHalfFiberLength = maxHalfFiberLength; - s.labels = labels; + s.seedLabels = seedLabels; + s.stopLabels = stopLabels; s.num_threads = numThreads; s.Qm = l_Qm; @@ -307,6 +311,10 @@ int ukf_parse_cli(int argc, char** argv, UKFSettings& s) s.output_file_with_second_tensor = tractsWithSecondTensor; s.dwiFile = dwiFile; s.seedsFile = seedsFile; + s.stopFile = stopFile; + s.wmFile = wmFile; + s.gmFile = gmFile; + s.csfFile = csfFile; s.maskFile = maskFile; s.writeAsciiTracts = writeAsciiTracts; s.writeUncompressedTracts = writeUncompressedTracts; diff --git a/ukf/tractography.cc b/ukf/tractography.cc index 9332019..e4e3c95 100644 --- a/ukf/tractography.cc +++ b/ukf/tractography.cc @@ -82,7 +82,11 @@ Tractography::Tractography(UKFSettings s) : _max_length(static_cast(std::ceil(s.maxHalfFiberLength / s.stepLength) ) ), _full_brain(false), _noddi(s.noddi), - _fa_min(s.fa_min), _mean_signal_min(s.mean_signal_min), + _wm_prob_threshold(s.wm_prob_threshold), + _gm_prob_threshold(s.gm_prob_threshold), + _csf_prob_threshold(s.csf_prob_threshold), + _fa_min(s.fa_min), + _mean_signal_min(s.mean_signal_min), _seeding_threshold(s.seeding_threshold), _num_tensors(s.num_tensors), _seeds_per_voxel(s.seeds_per_voxel), @@ -92,7 +96,8 @@ Tractography::Tractography(UKFSettings s) : _free_water(s.free_water), _stepLength(s.stepLength), _steps_per_record(s.recordLength/s.stepLength), - _labels(s.labels), + _seed_labels(s.seedLabels), + _stop_labels(s.stopLabels), Qm(s.Qm), Ql(s.Ql), @@ -316,8 +321,10 @@ void Tractography::UpdateFilterModelType() _model->set_signal_dim(_signal_data->GetSignalDimension() * 2); } -bool Tractography::SetData(void* data, void* mask, void* seed, - bool normalizedDWIData) +bool Tractography::SetData(void* data, void* mask, + bool normalizedDWIData, + void* seed, void* stop, + void* wm, void* gm, void* csf) { if (!data || !mask) { @@ -325,32 +332,35 @@ bool Tractography::SetData(void* data, void* mask, void* seed, return true; } - if (!seed) + if (!seed && !wm) // if no seeding mask or wm map is provided, brain mask with thresholds will be used for seeding { _full_brain = true; } _signal_data = new NrrdData(_sigma_signal, _sigma_mask); - _signal_data->SetData((Nrrd*)data, (Nrrd*)mask, (Nrrd*)seed, normalizedDWIData); + _signal_data->SetData((Nrrd*)data, (Nrrd*)mask, normalizedDWIData, (Nrrd*)seed, (Nrrd*)stop, (Nrrd*)wm, (Nrrd*)gm, (Nrrd*)csf); return false; } bool Tractography::LoadFiles(const std::string& data_file, - const std::string& seed_file, const std::string& mask_file, const bool normalized_DWI_data, - const bool output_normalized_DWI_data - ) + const bool output_normalized_DWI_data, + const std::string& seed_file, + const std::string& stop_file, + const std::string& wm_file, + const std::string& gm_file, + const std::string& csf_file) { _signal_data = new NrrdData(_sigma_signal, _sigma_mask); - if( seed_file.empty() ) + if( seed_file.empty() && wm_file.empty() ) { _full_brain = true; } - if( _signal_data->LoadData(data_file, seed_file, mask_file, normalized_DWI_data, output_normalized_DWI_data) ) + if( _signal_data->LoadData(data_file, mask_file, normalized_DWI_data, output_normalized_DWI_data, seed_file, stop_file, wm_file, gm_file, csf_file) ) { std::cout << "ISignalData could not be loaded" << std::endl; delete _signal_data; @@ -370,7 +380,7 @@ void Tractography::Init(std::vector& seed_infos) int signal_dim = _signal_data->GetSignalDimension(); stdVec_t seeds; - if (! (_labels.size() > 0)) + if (! (_seed_labels.size() > 0)) { std::cout << "No label data!"; throw; @@ -382,11 +392,12 @@ void Tractography::Init(std::vector& seed_infos) } else if(!_full_brain) { - _signal_data->GetSeeds(_labels, seeds); + _signal_data->GetSeeds(_seed_labels, _wm_prob_threshold, seeds); } else { // Iterate through all brain voxels and take those as seeds voxels. + std::cout << "Seeding using Option 1 (or no option has been specified), from the input brain mask." << std::endl; const vec3_t dim = _signal_data->dim(); for( int x = 0; x < dim[0]; ++x ) { @@ -1336,6 +1347,8 @@ void Tractography::Follow2T(const int thread_id, // the fiber gets too long. const bool is_brain = _signal_data->ScalarMaskValue(x) > 0; // _signal_data->Interp3ScalarMask(x) > 0.1; // is this 0.1 correct? yes + const bool is_stop = _signal_data->ScalarStopValue(_stop_labels, _gm_prob_threshold, _csf_prob_threshold, x) > 0; + // after here state does not change until next step. state_tmp.col(0) = state; @@ -1343,13 +1356,25 @@ void Tractography::Follow2T(const int thread_id, _model->H(state_tmp, signal_tmp); // signal_tmp is written, but only used to calculate mean signal const ukfPrecisionType mean_signal = s2adc(signal_tmp); - const bool in_csf = (_noddi) ? ( mean_signal < _mean_signal_min ) : - ( mean_signal < _mean_signal_min || fa < _fa_min); + + bool in_csf; + if (_signal_data->isGMCSFProvided()) + { + in_csf = false; + } + else + { + if(_noddi) + in_csf = mean_signal < _mean_signal_min; + else + in_csf = mean_signal < _mean_signal_min || fa < _fa_min; + } const bool is_curving = curve_radius(fiber.position) < _min_radius; if( !is_brain || in_csf + || is_stop || stepnr > _max_length // Stop when the fiber is too long || is_curving) { @@ -1463,6 +1488,9 @@ void Tractography::Follow1T(const int thread_id, // Terminate if off brain or in CSF. const bool is_brain = _signal_data->ScalarMaskValue(x) > 0; //_signal_data->Interp3ScalarMask(x) > 0.1; // x is the seed point + + const bool is_stop = _signal_data->ScalarStopValue(_stop_labels, _gm_prob_threshold, _csf_prob_threshold, x) > 0; + state_tmp.col(0) = state; _model->H(state_tmp, signal_tmp); @@ -1470,15 +1498,23 @@ void Tractography::Follow1T(const int thread_id, // Check mean_signal threshold const ukfPrecisionType mean_signal = s2adc(signal_tmp); bool in_csf; - if(_noddi) - in_csf = mean_signal < _mean_signal_min; + if (_signal_data->isGMCSFProvided()) + { + in_csf = false; + } else - in_csf = mean_signal < _mean_signal_min || fa < _fa_min; + { + if(_noddi) + in_csf = mean_signal < _mean_signal_min; + else + in_csf = mean_signal < _mean_signal_min || fa < _fa_min; + } bool is_curving = curve_radius(fiber.position) < _min_radius; if( !is_brain || in_csf + || is_stop || stepnr > _max_length // Stop when fiber is too long || is_curving ) { diff --git a/ukf/tractography.h b/ukf/tractography.h index c75fab3..4a4fe2e 100644 --- a/ukf/tractography.h +++ b/ukf/tractography.h @@ -38,6 +38,9 @@ struct UKFSettings { bool transform_position; bool store_glyphs; bool branches_only; + ukfPrecisionType wm_prob_threshold; + ukfPrecisionType gm_prob_threshold; + ukfPrecisionType csf_prob_threshold; ukfPrecisionType fa_min; ukfPrecisionType mean_signal_min; ukfPrecisionType seeding_threshold; @@ -51,7 +54,8 @@ struct UKFSettings { ukfPrecisionType stepLength; ukfPrecisionType recordLength; ukfPrecisionType maxHalfFiberLength; - std::vector labels; + std::vector seedLabels; + std::vector stopLabels; ukfPrecisionType Qm; ukfPrecisionType Ql; @@ -78,6 +82,10 @@ struct UKFSettings { std::string output_file_with_second_tensor; std::string dwiFile; std::string seedsFile; + std::string stopFile; + std::string wmFile; + std::string gmFile; + std::string csfFile; std::string maskFile; }; @@ -115,14 +123,16 @@ friend class vtkSlicerInteractiveUKFLogic; * Load the files that contain the DWI signal, the seeds and a mask * defining the volume of the brain. */ - bool LoadFiles(const std::string& data_file, const std::string& seed_file, const std::string& mask_file, - const bool normalized_DWI_data, const bool output_normalized_DWI_data); + bool LoadFiles(const std::string& data_file, const std::string& mask_file, + const bool normalized_DWI_data, const bool output_normalized_DWI_data, + const std::string& seed_file, const std::string& stop_file, + const std::string& wm_file, const std::string& gm_file, const std::string& csf_file); /** * Directly set the data volume pointers */ - bool SetData(void* data, void* mask, void* seed, bool normalizedDWIData); + bool SetData(void* data, void* mask, bool normalizedDWIData, void* seed, void* stop, void* wm, void* gm, void* csf); /** * Directly set the seed locations @@ -294,6 +304,9 @@ friend class vtkSlicerInteractiveUKFLogic; int _nPosFreeWater; // Parameters for the tractography + ukfPrecisionType _wm_prob_threshold; + ukfPrecisionType _gm_prob_threshold; + ukfPrecisionType _csf_prob_threshold; ukfPrecisionType _fa_min; ukfPrecisionType _mean_signal_min; ukfPrecisionType _seeding_threshold; @@ -305,7 +318,8 @@ friend class vtkSlicerInteractiveUKFLogic; bool _free_water; ukfPrecisionType _stepLength; int _steps_per_record; - std::vector _labels; + std::vector _seed_labels; + std::vector _stop_labels; stdVec_t _ext_seeds; ukfPrecisionType Qm;