From 9248b5b05005cf15d52dbf980985a35a18bbd589 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 27 Jun 2017 15:17:07 +1000 Subject: [PATCH 01/11] fixelcfestats: New option -mask This option restricts the statistical analysis to a subset of fixels, as defined by a boolean fixel data file. Calculation of the fixel-fixel connectivity matrix, data smoothing, t-statistic calculation, and stiatistical enhancement, are all only performed on those fixels within the mask. However, the output fixel directory contains all of the fixels from the input template fixel image, not just those within the mask; this retains fixel-fixel correspondence between the input template image and the output results. For fixels outside the mask, a value of NaN is written for all fixel data file outputs. --- cmd/fixelcfestats.cpp | 205 +++++++++++++++++++++++++++--------------- src/stats/cfe.cpp | 18 ++-- src/stats/cfe.h | 4 +- 3 files changed, 149 insertions(+), 78 deletions(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index b6e60bfbe7..5813d18a09 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -109,7 +109,10 @@ void usage () + Argument ("threshold").type_float (0.0, 1.0) + Option ("angle", "the max angle threshold for assigning streamline tangents to fixels (Default: " + str(DEFAULT_ANGLE_THRESHOLD, 2) + " degrees)") - + Argument ("value").type_float (0.0, 90.0); + + Argument ("value").type_float (0.0, 90.0) + + + Option ("mask", "provide a fixel data file containing a mask of those fixels to be used during processing") + + Argument ("file").type_image_in(); } @@ -118,12 +121,14 @@ void usage () template void write_fixel_output (const std::string& filename, const VectorType& data, + const vector& fixel2row, const Header& header) { + assert (header.size(0) == fixel2row.size()); auto output = Image::create (filename, header); - for (uint32_t i = 0; i < data.size(); ++i) { + for (uint32_t i = 0; i != fixel2row.size(); ++i) { output.index(0) = i; - output.value() = data[i]; + output.value() = (fixel2row[i] >= 0) ? data[fixel2row[i]] : NaN; } } @@ -148,10 +153,37 @@ void run() { const std::string input_fixel_directory = argument[0]; Header index_header = Fixel::find_index_header (input_fixel_directory); auto index_image = index_header.get_image(); - const uint32_t num_fixels = Fixel::get_number_of_fixels (index_header); CONSOLE ("number of fixels: " + str(num_fixels)); + Image mask; + uint32_t mask_fixels; + // Lookup tables that will map from input fixel index to row number, and + // from row number to fixel index + vector fixel2row (num_fixels); + opt = get_options ("mask"); + if (opt.size()) { + mask = Image::open (opt[0][0]); + Fixel::check_data_file (mask); + if (!Fixel::fixels_match (index_header, mask)) + throw Exception ("Mask image provided using -mask option does not match fixel template"); + uint32_t mask_fixels = 0; + for (auto l = Loop(mask) (mask); l; ++l) + fixel2row[mask.index(0)] = mask.value() ? mask_fixels++ : -1; + CONSOLE ("Template mask contains " + str(mask_fixels) + " fixels"); + } else { + Header data_header; + data_header.ndim() = 1; + data_header.size(0) = num_fixels; + data_header.spacing(0) = NaN; + mask = Image::scratch (data_header, "scratch fixel mask"); + for (auto l = Loop(mask) (mask); l; ++l) + mask.value() = true; + mask_fixels = num_fixels; + for (uint32_t f = 0; f != num_fixels; ++f) + fixel2row[f] = f; + } + vector positions (num_fixels); vector directions (num_fixels); @@ -163,13 +195,14 @@ void run() { // Load template fixel directions Transform image_transform (index_image); for (auto i = Loop ("loading template fixel directions and positions", index_image, 0, 3)(index_image); i; ++i) { - const Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2)); + Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2)); + vox = image_transform.voxel2scanner * vox; index_image.index(3) = 1; uint32_t offset = index_image.value(); size_t fixel_index = 0; for (auto f = Fixel::Loop (index_image) (directions_data); f; ++f, ++fixel_index) { directions[offset + fixel_index] = directions_data.row(1); - positions[offset + fixel_index] = image_transform.voxel2scanner * vox; + positions[offset + fixel_index] = vox; } } } @@ -229,36 +262,34 @@ void run() { // Compute fixel-fixel connectivity vector > connectivity_matrix (num_fixels); - vector fixel_TDI (num_fixels, 0.0); + vector fixel_TDI (num_fixels, 0); const std::string track_filename = argument[4]; DWI::Tractography::Properties properties; DWI::Tractography::Reader track_file (track_filename, properties); // Read in tracts, and compute whole-brain fixel-fixel connectivity - const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); + const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); if (!num_tracks) throw Exception ("no tracks found in input file"); if (num_tracks < 1000000) WARN ("more than 1 million tracks should be used to ensure robust fixel-fixel connectivity"); { - using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir; DWI::Tractography::Mapping::TrackLoader loader (track_file, num_tracks, "pre-computing fixel-fixel connectivity"); DWI::Tractography::Mapping::TrackMapperBase mapper (index_image); mapper.set_upsample_ratio (DWI::Tractography::Mapping::determine_upsample_ratio (index_header, properties, 0.333f)); mapper.set_use_precise_mapping (true); - Stats::CFE::TrackProcessor tract_processor (index_image, directions, fixel_TDI, connectivity_matrix, angular_threshold); + Stats::CFE::TrackProcessor tract_processor (index_image, directions, mask, fixel_TDI, connectivity_matrix, angular_threshold); Thread::run_queue ( loader, Thread::batch (DWI::Tractography::Streamline()), mapper, - Thread::batch (SetVoxelDir()), + Thread::batch (DWI::Tractography::Mapping::SetVoxelDir()), tract_processor); } track_file.close(); // Normalise connectivity matrix and threshold, pre-compute fixel-fixel weights for smoothing. - vector > smoothing_weights (num_fixels); + vector > smoothing_weights (mask_fixels); bool do_smoothing = false; - const float gaussian_const2 = 2.0 * smooth_std_dev * smooth_std_dev; float gaussian_const1 = 1.0; if (smooth_std_dev > 0.0) { @@ -269,43 +300,71 @@ void run() { { ProgressBar progress ("normalising and thresholding fixel-fixel connectivity matrix", num_fixels); for (uint32_t fixel = 0; fixel < num_fixels; ++fixel) { - - auto it = connectivity_matrix[fixel].begin(); - while (it != connectivity_matrix[fixel].end()) { - const connectivity_value_type connectivity = it->second.value / connectivity_value_type (fixel_TDI[fixel]); - if (connectivity < connectivity_threshold) { - connectivity_matrix[fixel].erase (it++); - } else { - if (do_smoothing) { - const value_type distance = std::sqrt (Math::pow2 (positions[fixel][0] - positions[it->first][0]) + - Math::pow2 (positions[fixel][1] - positions[it->first][1]) + - Math::pow2 (positions[fixel][2] - positions[it->first][2])); - const connectivity_value_type smoothing_weight = connectivity * gaussian_const1 * std::exp (-std::pow (distance, 2) / gaussian_const2); - if (smoothing_weight > 0.01) - smoothing_weights[fixel].insert (std::pair (it->first, smoothing_weight)); + mask.index(0) = fixel; + const uint32_t row = fixel2row[fixel]; + if (mask.value()) { + + // Here, the connectivity matrix needs to be modified to reflect the + // fact that fixel indices in the template fixel image may not + // correspond to rows in the statistical analysis + std::map new_connectivity_matrix_row; + std::map smoothing_weights_row; + connectivity_value_type sum_weights = 0.0; + + for (auto& it : connectivity_matrix[fixel]) { +#ifndef NDEBUG + // Even if this fixel is within the mask, it should still not + // connect to any fixel that is outside the mask + mask.index(0) = it.first; + assert (!mask.value()); +#endif + const connectivity_value_type connectivity = it.second.value / connectivity_value_type (fixel_TDI[fixel]); + if (connectivity >= connectivity_threshold) { + if (do_smoothing) { + const value_type distance = std::sqrt (Math::pow2 (positions[fixel][0] - positions[it.first][0]) + + Math::pow2 (positions[fixel][1] - positions[it.first][1]) + + Math::pow2 (positions[fixel][2] - positions[it.first][2])); + const connectivity_value_type smoothing_weight = connectivity * gaussian_const1 * std::exp (-std::pow (distance, 2) / gaussian_const2); + if (smoothing_weight >= connectivity_threshold) { + smoothing_weights_row[fixel2row[it.first]] = smoothing_weight; + sum_weights += smoothing_weight; + } + } + // Here we pre-exponentiate each connectivity value by C + new_connectivity_matrix_row[fixel2row[it.first]] = std::pow (connectivity, cfe_c); } - // Here we pre-exponentiate each connectivity value by C - it->second.value = std::pow (connectivity, cfe_c); - ++it; } - } - // Make sure the fixel is fully connected to itself - connectivity_matrix[fixel].insert (std::pair (fixel, Stats::CFE::connectivity (1.0))); - smoothing_weights[fixel].insert (std::pair (fixel, gaussian_const1)); - - // Normalise smoothing weights - value_type sum = 0.0; - for (auto smooth_it = smoothing_weights[fixel].begin(); smooth_it != smoothing_weights[fixel].end(); ++smooth_it) { - sum += smooth_it->second; - } - value_type norm_factor = 1.0 / sum; - for (auto smooth_it = smoothing_weights[fixel].begin(); smooth_it != smoothing_weights[fixel].end(); ++smooth_it) { - smooth_it->second *= norm_factor; + + // Make sure the fixel is fully connected to itself + new_connectivity_matrix_row[row] = Stats::CFE::connectivity (1.0); + smoothing_weights_row[row] = connectivity_value_type (gaussian_const1); + sum_weights += gaussian_const1; + + // Normalise smoothing weights + const connectivity_value_type norm_factor = 1.0 / sum_weights; + for (auto i : smoothing_weights_row) + i.second *= norm_factor; + + // Write the new data for this fixel to the relevant structures + // Safe to over-write connectivity_matrix, since the row in the + // statistical analysis corresponding to this fixel is always + // less than or equal to the fixel index. + connectivity_matrix[row] = std::move (new_connectivity_matrix_row); + smoothing_weights[row] = std::move (smoothing_weights_row); + + } else { + // If fixel is not in the mask, tract_processor should never assign + // any connections to it + assert (connectivity_matrix[fixel].empty()); } progress++; } } + // If the connectivity matrix has shrunk in size due to the masking of fixels, + // throw out any now-erroneous data off the end of the vector + connectivity_matrix.resize (mask_fixels); + Header output_header (header); output_header.keyval()["num permutations"] = str(num_perms); output_header.keyval()["dh"] = str(cfe_dh); @@ -318,7 +377,7 @@ void run() { // Load input data - matrix_type data (num_fixels, identifiers.size()); + matrix_type data (mask_fixels, identifiers.size()); data.setZero(); { ProgressBar progress ("loading input images", identifiers.size()); @@ -326,7 +385,7 @@ void run() { LogLevelLatch log_level (0); auto subject_data = Image::open (identifiers[subject]).with_direct_io(); - vector subject_data_vector (num_fixels, 0.0); + vector subject_data_vector (mask_fixels, 0.0); for (auto i = Loop (index_image, 0, 3)(index_image); i; ++i) { index_image.index(3) = 1; uint32_t offset = index_image.value(); @@ -334,18 +393,22 @@ void run() { for (auto f = Fixel::Loop (index_image) (subject_data); f; ++f, ++fixel_index) { if (!std::isfinite(static_cast(subject_data.value()))) throw Exception ("subject data file " + identifiers[subject] + " contains non-finite value: " + str(subject_data.value())); - subject_data_vector[offset + fixel_index] = subject_data.value(); + // Note: Data mapped from fixel index to row here; + // smoothing weights are prepared in this format + subject_data_vector[fixel2row[offset+fixel_index]] = subject_data.value(); } } // Smooth the data - for (size_t fixel = 0; fixel < num_fixels; ++fixel) { - value_type value = 0.0; - std::map::const_iterator it = smoothing_weights[fixel].begin(); - for (; it != smoothing_weights[fixel].end(); ++it) { - value += subject_data_vector[it->first] * it->second; + if (do_smoothing) { + for (size_t fixel = 0; fixel < mask_fixels; ++fixel) { + value_type value = 0.0; + for (auto i : smoothing_weights[fixel]) + value += subject_data_vector[i.first] * i.second; + data (fixel, subject) = value; } - data (fixel, subject) = value; + } else { + data.col (subject) = subject_data_vector; } progress++; } @@ -359,15 +422,15 @@ void run() { auto temp = Math::Stats::GLM::solve_betas (data, design); for (ssize_t i = 0; i < contrast.cols(); ++i) { - write_fixel_output (Path::join (output_fixel_directory, "beta" + str(i) + ".mif"), temp.row(i), output_header); + write_fixel_output (Path::join (output_fixel_directory, "beta" + str(i) + ".mif"), temp.row(i), fixel2row, output_header); ++progress; } temp = Math::Stats::GLM::abs_effect_size (data, design, contrast); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "abs_effect.mif"), temp.row(0), output_header); ++progress; + write_fixel_output (Path::join (output_fixel_directory, "abs_effect.mif"), temp.row(0), fixel2row, output_header); ++progress; temp = Math::Stats::GLM::std_effect_size (data, design, contrast); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "std_effect.mif"), temp.row(0), output_header); ++progress; + write_fixel_output (Path::join (output_fixel_directory, "std_effect.mif"), temp.row(0), fixel2row, output_header); ++progress; temp = Math::Stats::GLM::stdev (data, design); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "std_dev.mif"), temp.row(0), output_header); + write_fixel_output (Path::join (output_fixel_directory, "std_dev.mif"), temp.row(0), fixel2row, output_header); } Math::Stats::GLMTTest glm_ttest (data, design, contrast); @@ -386,35 +449,35 @@ void run() { Stats::PermTest::precompute_empirical_stat (glm_ttest, cfe_integrator, permutations, empirical_cfe_statistic); } output_header.keyval()["nonstationary adjustment"] = str(true); - write_fixel_output (Path::join (output_fixel_directory, "cfe_empirical.mif"), empirical_cfe_statistic, output_header); + write_fixel_output (Path::join (output_fixel_directory, "cfe_empirical.mif"), empirical_cfe_statistic, fixel2row, output_header); } else { output_header.keyval()["nonstationary adjustment"] = str(false); } // Precompute default statistic and CFE statistic - vector_type cfe_output (num_fixels); + vector_type cfe_output (mask_fixels); std::shared_ptr cfe_output_neg; - vector_type tvalue_output (num_fixels); + vector_type tvalue_output (mask_fixels); if (compute_negative_contrast) - cfe_output_neg.reset (new vector_type (num_fixels)); + cfe_output_neg.reset (new vector_type (mask_fixels)); Stats::PermTest::precompute_default_permutation (glm_ttest, cfe_integrator, empirical_cfe_statistic, cfe_output, cfe_output_neg, tvalue_output); - write_fixel_output (Path::join (output_fixel_directory, "cfe.mif"), cfe_output, output_header); - write_fixel_output (Path::join (output_fixel_directory, "tvalue.mif"), tvalue_output, output_header); + write_fixel_output (Path::join (output_fixel_directory, "cfe.mif"), cfe_output, fixel2row, output_header); + write_fixel_output (Path::join (output_fixel_directory, "tvalue.mif"), tvalue_output, fixel2row, output_header); if (compute_negative_contrast) - write_fixel_output (Path::join (output_fixel_directory, "cfe_neg.mif"), *cfe_output_neg, output_header); + write_fixel_output (Path::join (output_fixel_directory, "cfe_neg.mif"), *cfe_output_neg, fixel2row, output_header); // Perform permutation testing if (!get_options ("notest").size()) { vector_type perm_distribution (num_perms); std::shared_ptr perm_distribution_neg; - vector_type uncorrected_pvalues (num_fixels); + vector_type uncorrected_pvalues (mask_fixels); std::shared_ptr uncorrected_pvalues_neg; if (compute_negative_contrast) { perm_distribution_neg.reset (new vector_type (num_perms)); - uncorrected_pvalues_neg.reset (new vector_type (num_fixels)); + uncorrected_pvalues_neg.reset (new vector_type (mask_fixels)); } if (permutations.size()) { @@ -432,17 +495,17 @@ void run() { ProgressBar progress ("outputting final results"); save_matrix (perm_distribution, Path::join (output_fixel_directory, "perm_dist.txt")); ++progress; - vector_type pvalue_output (num_fixels); + vector_type pvalue_output (mask_fixels); Math::Stats::Permutation::statistic2pvalue (perm_distribution, cfe_output, pvalue_output); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "fwe_pvalue.mif"), pvalue_output, output_header); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "uncorrected_pvalue.mif"), uncorrected_pvalues, output_header); ++progress; + write_fixel_output (Path::join (output_fixel_directory, "fwe_pvalue.mif"), pvalue_output, fixel2row, output_header); ++progress; + write_fixel_output (Path::join (output_fixel_directory, "uncorrected_pvalue.mif"), uncorrected_pvalues, fixel2row, output_header); ++progress; if (compute_negative_contrast) { save_matrix (*perm_distribution_neg, Path::join (output_fixel_directory, "perm_dist_neg.txt")); ++progress; - vector_type pvalue_output_neg (num_fixels); + vector_type pvalue_output_neg (mask_fixels); Math::Stats::Permutation::statistic2pvalue (*perm_distribution_neg, *cfe_output_neg, pvalue_output_neg); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "fwe_pvalue_neg.mif"), pvalue_output_neg, output_header); ++progress; - write_fixel_output (Path::join (output_fixel_directory, "uncorrected_pvalue_neg.mif"), *uncorrected_pvalues_neg, output_header); + write_fixel_output (Path::join (output_fixel_directory, "fwe_pvalue_neg.mif"), pvalue_output_neg, fixel2row, output_header); ++progress; + write_fixel_output (Path::join (output_fixel_directory, "uncorrected_pvalue_neg.mif"), *uncorrected_pvalues_neg, fixel2row, output_header); } } } diff --git a/src/stats/cfe.cpp b/src/stats/cfe.cpp index 79c6bb017b..1f012fda53 100644 --- a/src/stats/cfe.cpp +++ b/src/stats/cfe.cpp @@ -25,11 +25,13 @@ namespace MR TrackProcessor::TrackProcessor (Image& fixel_indexer, const vector& fixel_directions, + Image& fixel_mask, vector& fixel_TDI, vector >& connectivity_matrix, const value_type angular_threshold) : fixel_indexer (fixel_indexer) , fixel_directions (fixel_directions), + fixel_mask (fixel_mask), fixel_TDI (fixel_TDI), connectivity_matrix (connectivity_matrix), angular_threshold_dp (std::cos (angular_threshold * (Math::pi/180.0))) { } @@ -43,22 +45,26 @@ namespace MR for (SetVoxelDir::const_iterator i = in.begin(); i != in.end(); ++i) { assign_pos_of (*i).to (fixel_indexer); fixel_indexer.index(3) = 0; - uint32_t num_fibres = fixel_indexer.value(); - if (num_fibres > 0) { + uint32_t num_fixels = fixel_indexer.value(); + if (num_fixels > 0) { fixel_indexer.index(3) = 1; uint32_t first_index = fixel_indexer.value(); - uint32_t last_index = first_index + num_fibres; - uint32_t closest_fixel_index = 0; + uint32_t last_index = first_index + num_fixels; + // Note: Streamlines can still be assigned to a fixel that is outside the mask; + // however this will not be permitted to contribute to the matrix + uint32_t closest_fixel_index = num_fixels; value_type largest_dp = 0.0; const direction_type dir (i->get_dir().normalized()); for (uint32_t j = first_index; j < last_index; ++j) { const value_type dp = std::abs (dir.dot (fixel_directions[j])); if (dp > largest_dp) { largest_dp = dp; - closest_fixel_index = j; + fixel_mask.index(0) = j; + if (fixel_mask.value()) + closest_fixel_index = j; } } - if (largest_dp > angular_threshold_dp) { + if (closest_fixel_index != num_fixels && largest_dp > angular_threshold_dp) { tract_fixel_indices.push_back (closest_fixel_index); fixel_TDI[closest_fixel_index]++; } diff --git a/src/stats/cfe.h b/src/stats/cfe.h index 9560e30756..1130be419e 100644 --- a/src/stats/cfe.h +++ b/src/stats/cfe.h @@ -43,7 +43,7 @@ namespace MR @{ */ - class connectivity { MEMALIGN(connectivity) + class connectivity { NOMEMALIGN public: connectivity () : value (0.0) { } connectivity (const connectivity_value_type v) : value (v) { } @@ -61,6 +61,7 @@ namespace MR public: TrackProcessor (Image& fixel_indexer, const vector& fixel_directions, + Image& fixel_mask, vector& fixel_TDI, vector >& connectivity_matrix, const value_type angular_threshold); @@ -70,6 +71,7 @@ namespace MR private: Image fixel_indexer; const vector& fixel_directions; + Image fixel_mask; vector& fixel_TDI; vector >& connectivity_matrix; const value_type angular_threshold_dp; From 11911d5ff34bf8ed54c31b59d30f325a5435e94f Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 27 Jun 2017 15:36:51 +1000 Subject: [PATCH 02/11] Docs update --- docs/reference/commands/fixelcfestats.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/reference/commands/fixelcfestats.rst b/docs/reference/commands/fixelcfestats.rst index f04ed4dfdb..209c92d544 100644 --- a/docs/reference/commands/fixelcfestats.rst +++ b/docs/reference/commands/fixelcfestats.rst @@ -62,6 +62,8 @@ Additional options for fixelcfestats - **-angle value** the max angle threshold for assigning streamline tangents to fixels (Default: 45 degrees) +- **-mask file** provide a fixel data file containing a mask of those fixels to be used during processing + Standard options ^^^^^^^^^^^^^^^^ From 7c58e1c131cf6d942279792cf2fbf68274fe70bd Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 27 Jun 2017 15:45:37 +1000 Subject: [PATCH 03/11] fixelcfestats: Fix compilation --- cmd/fixelcfestats.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 5813d18a09..54a47e3be8 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -408,7 +408,7 @@ void run() { data (fixel, subject) = value; } } else { - data.col (subject) = subject_data_vector; + data.col (subject) = Eigen::Map > (subject_data_vector.data(), mask_fixels); } progress++; } From 3a6aded2603916acabecc9eb5e5eb8d0d3389a0e Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 28 Jun 2017 15:12:34 +1000 Subject: [PATCH 04/11] fixelcfestats: Various fixels for -mask option Also includes a bug fix where opting to not smooth the input data would likely have resulted in an error. --- cmd/fixelcfestats.cpp | 59 +++++++++++++++++++++++++------------------ 1 file changed, 35 insertions(+), 24 deletions(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 54a47e3be8..d446cf83a9 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -124,7 +124,7 @@ void write_fixel_output (const std::string& filename, const vector& fixel2row, const Header& header) { - assert (header.size(0) == fixel2row.size()); + assert (size_t(header.size(0)) == fixel2row.size()); auto output = Image::create (filename, header); for (uint32_t i = 0; i != fixel2row.size(); ++i) { output.index(0) = i; @@ -153,13 +153,15 @@ void run() { const std::string input_fixel_directory = argument[0]; Header index_header = Fixel::find_index_header (input_fixel_directory); auto index_image = index_header.get_image(); + const uint32_t num_fixels = Fixel::get_number_of_fixels (index_header); CONSOLE ("number of fixels: " + str(num_fixels)); Image mask; uint32_t mask_fixels; - // Lookup tables that will map from input fixel index to row number, and - // from row number to fixel index + // Lookup table that maps from input fixel index to row number + // Note that if a fixel is masked out, it will have a value of -1 + // in this array vector fixel2row (num_fixels); opt = get_options ("mask"); if (opt.size()) { @@ -170,15 +172,20 @@ void run() { uint32_t mask_fixels = 0; for (auto l = Loop(mask) (mask); l; ++l) fixel2row[mask.index(0)] = mask.value() ? mask_fixels++ : -1; - CONSOLE ("Template mask contains " + str(mask_fixels) + " fixels"); + CONSOLE ("Mask contains " + str(mask_fixels) + " fixels"); } else { Header data_header; - data_header.ndim() = 1; + data_header.ndim() = 3; data_header.size(0) = num_fixels; - data_header.spacing(0) = NaN; + data_header.size(1) = 1; + data_header.size(2) = 1; + data_header.spacing(0) = data_header.spacing(1) = data_header.spacing(2) = NaN; + data_header.stride(0) = 1; data_header.stride(1) = 2; data_header.stride(2) = 3; mask = Image::scratch (data_header, "scratch fixel mask"); - for (auto l = Loop(mask) (mask); l; ++l) + for (uint32_t f = 0; f != num_fixels; ++f) { + mask.index(0) = f; mask.value() = true; + } mask_fixels = num_fixels; for (uint32_t f = 0; f != num_fixels; ++f) fixel2row[f] = f; @@ -290,6 +297,7 @@ void run() { // Normalise connectivity matrix and threshold, pre-compute fixel-fixel weights for smoothing. vector > smoothing_weights (mask_fixels); bool do_smoothing = false; + const float gaussian_const2 = 2.0 * smooth_std_dev * smooth_std_dev; float gaussian_const1 = 1.0; if (smooth_std_dev > 0.0) { @@ -301,14 +309,14 @@ void run() { ProgressBar progress ("normalising and thresholding fixel-fixel connectivity matrix", num_fixels); for (uint32_t fixel = 0; fixel < num_fixels; ++fixel) { mask.index(0) = fixel; - const uint32_t row = fixel2row[fixel]; + const int32_t row = fixel2row[fixel]; + if (mask.value()) { // Here, the connectivity matrix needs to be modified to reflect the // fact that fixel indices in the template fixel image may not // correspond to rows in the statistical analysis std::map new_connectivity_matrix_row; - std::map smoothing_weights_row; connectivity_value_type sum_weights = 0.0; for (auto& it : connectivity_matrix[fixel]) { @@ -316,7 +324,7 @@ void run() { // Even if this fixel is within the mask, it should still not // connect to any fixel that is outside the mask mask.index(0) = it.first; - assert (!mask.value()); + assert (mask.value()); #endif const connectivity_value_type connectivity = it.second.value / connectivity_value_type (fixel_TDI[fixel]); if (connectivity >= connectivity_threshold) { @@ -326,7 +334,7 @@ void run() { Math::pow2 (positions[fixel][2] - positions[it.first][2])); const connectivity_value_type smoothing_weight = connectivity * gaussian_const1 * std::exp (-std::pow (distance, 2) / gaussian_const2); if (smoothing_weight >= connectivity_threshold) { - smoothing_weights_row[fixel2row[it.first]] = smoothing_weight; + smoothing_weights[row][fixel2row[it.first]] = smoothing_weight; sum_weights += smoothing_weight; } } @@ -337,26 +345,27 @@ void run() { // Make sure the fixel is fully connected to itself new_connectivity_matrix_row[row] = Stats::CFE::connectivity (1.0); - smoothing_weights_row[row] = connectivity_value_type (gaussian_const1); - sum_weights += gaussian_const1; + smoothing_weights[row][row] = connectivity_value_type(gaussian_const1); + sum_weights += connectivity_value_type(gaussian_const1); // Normalise smoothing weights - const connectivity_value_type norm_factor = 1.0 / sum_weights; - for (auto i : smoothing_weights_row) + const connectivity_value_type norm_factor = connectivity_value_type(1.0) / sum_weights; + for (auto i : smoothing_weights[row]) i.second *= norm_factor; - // Write the new data for this fixel to the relevant structures - // Safe to over-write connectivity_matrix, since the row in the - // statistical analysis corresponding to this fixel is always - // less than or equal to the fixel index. + // Write the new connectivity matrix data to the relevant structure + // Note that due to fixel masking, this may not be the same row as + // the fixel index connectivity_matrix[row] = std::move (new_connectivity_matrix_row); - smoothing_weights[row] = std::move (smoothing_weights_row); - + // Force deallocation of memory used + std::map().swap (new_connectivity_matrix_row); } else { // If fixel is not in the mask, tract_processor should never assign // any connections to it assert (connectivity_matrix[fixel].empty()); + } + progress++; } } @@ -365,6 +374,7 @@ void run() { // throw out any now-erroneous data off the end of the vector connectivity_matrix.resize (mask_fixels); + Header output_header (header); output_header.keyval()["num permutations"] = str(num_perms); output_header.keyval()["dh"] = str(cfe_dh); @@ -393,9 +403,10 @@ void run() { for (auto f = Fixel::Loop (index_image) (subject_data); f; ++f, ++fixel_index) { if (!std::isfinite(static_cast(subject_data.value()))) throw Exception ("subject data file " + identifiers[subject] + " contains non-finite value: " + str(subject_data.value())); - // Note: Data mapped from fixel index to row here; - // smoothing weights are prepared in this format - subject_data_vector[fixel2row[offset+fixel_index]] = subject_data.value(); + // Note that immediately on import, data are re-arranged according to fixel mask + const int32_t row = fixel2row[offset+fixel_index]; + if (row >= 0) + subject_data_vector[row] = subject_data.value(); } } From 36224e5445fbd88ed50c14a8c3d4b8f774e0273a Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Thu, 29 Jun 2017 11:32:30 +1000 Subject: [PATCH 05/11] fixelcfestats: Change type of normalised connectivity matrix On testing prior changes to fixelcfestats to support the -mask option, it was found to run extremely slowly, even when not using the -mask option. The suspected cause was that during the matrix normalisation, each row of the connectivity matrix (a map) was being re-constructed from scratch (due to the possibility of fixel indices changing), and hence what was a balanced binary search tree would become essentially a sorted doubly-linked-list. Now, each row of the connectivity matrix is instead converted to essentially a vector>, since it is no longer necessary to have fast lookup for an arbitrary fixel index. This should reduce memory requirement after matrix construction has completed, and will also hopefully improve cache performance since the connectivity data for each fixel will be serialised. There are also various little fixes and tweaks to polish off additions for supporting the -mask option. --- cmd/fixelcfestats.cpp | 76 ++++++++++++++++++++++--------------------- src/stats/cfe.cpp | 30 ++++++++--------- src/stats/cfe.h | 39 ++++++++++++++++++---- 3 files changed, 86 insertions(+), 59 deletions(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index d446cf83a9..ac06fb7ac6 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -39,6 +39,7 @@ using namespace MR::DWI::Tractography::Mapping; using namespace MR::Math::Stats; using Stats::CFE::direction_type; using Stats::CFE::connectivity_value_type; +using Stats::CFE::index_type; #define DEFAULT_CFE_DH 0.1 #define DEFAULT_CFE_E 2.0 @@ -50,7 +51,7 @@ using Stats::CFE::connectivity_value_type; void usage () { - AUTHOR = "David Raffelt (david.raffelt@florey.edu.au)"; + AUTHOR = "David Raffelt (david.raffelt@florey.edu.au) and Robert E. Smith (robert.smith@florey.edu.au)"; SYNOPSIS = "Fixel-based analysis using connectivity-based fixel enhancement and non-parametric permutation testing"; @@ -126,7 +127,7 @@ void write_fixel_output (const std::string& filename, { assert (size_t(header.size(0)) == fixel2row.size()); auto output = Image::create (filename, header); - for (uint32_t i = 0; i != fixel2row.size(); ++i) { + for (size_t i = 0; i != fixel2row.size(); ++i) { output.index(0) = i; output.value() = (fixel2row[i] >= 0) ? data[fixel2row[i]] : NaN; } @@ -152,16 +153,16 @@ void run() { const std::string input_fixel_directory = argument[0]; Header index_header = Fixel::find_index_header (input_fixel_directory); - auto index_image = index_header.get_image(); + auto index_image = index_header.get_image(); - const uint32_t num_fixels = Fixel::get_number_of_fixels (index_header); + const index_type num_fixels = Fixel::get_number_of_fixels (index_header); CONSOLE ("number of fixels: " + str(num_fixels)); Image mask; - uint32_t mask_fixels; + index_type mask_fixels; // Lookup table that maps from input fixel index to row number // Note that if a fixel is masked out, it will have a value of -1 - // in this array + // in this array; hence require a signed integer type vector fixel2row (num_fixels); opt = get_options ("mask"); if (opt.size()) { @@ -169,8 +170,8 @@ void run() { Fixel::check_data_file (mask); if (!Fixel::fixels_match (index_header, mask)) throw Exception ("Mask image provided using -mask option does not match fixel template"); - uint32_t mask_fixels = 0; - for (auto l = Loop(mask) (mask); l; ++l) + mask_fixels = 0; + for (mask.index(0) = 0; mask.index(0) != num_fixels; ++mask.index(0)) fixel2row[mask.index(0)] = mask.value() ? mask_fixels++ : -1; CONSOLE ("Mask contains " + str(mask_fixels) + " fixels"); } else { @@ -182,10 +183,8 @@ void run() { data_header.spacing(0) = data_header.spacing(1) = data_header.spacing(2) = NaN; data_header.stride(0) = 1; data_header.stride(1) = 2; data_header.stride(2) = 3; mask = Image::scratch (data_header, "scratch fixel mask"); - for (uint32_t f = 0; f != num_fixels; ++f) { - mask.index(0) = f; + for (mask.index(0) = 0; mask.index(0) != num_fixels; ++mask.index(0)) mask.value() = true; - } mask_fixels = num_fixels; for (uint32_t f = 0; f != num_fixels; ++f) fixel2row[f] = f; @@ -205,7 +204,7 @@ void run() { Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2)); vox = image_transform.voxel2scanner * vox; index_image.index(3) = 1; - uint32_t offset = index_image.value(); + const index_type offset = index_image.value(); size_t fixel_index = 0; for (auto f = Fixel::Loop (index_image) (directions_data); f; ++f, ++fixel_index) { directions[offset + fixel_index] = directions_data.row(1); @@ -217,7 +216,7 @@ void run() { vector identifiers; Header header; { - ProgressBar progress ("validating input files..."); + ProgressBar progress ("validating input files"); std::ifstream ifs (argument[1].c_str()); std::string temp; while (getline (ifs, temp)) { @@ -268,7 +267,7 @@ void run() { throw Exception ("only a single contrast vector (defined as a row) is currently supported"); // Compute fixel-fixel connectivity - vector > connectivity_matrix (num_fixels); + Stats::CFE::init_connectivity_matrix_type connectivity_matrix (num_fixels); vector fixel_TDI (num_fixels, 0); const std::string track_filename = argument[4]; DWI::Tractography::Properties properties; @@ -294,8 +293,10 @@ void run() { } track_file.close(); - // Normalise connectivity matrix and threshold, pre-compute fixel-fixel weights for smoothing. - vector > smoothing_weights (mask_fixels); + // Normalise connectivity matrix, threshold, and put in a more efficient format + Stats::CFE::norm_connectivity_matrix_type norm_connectivity_matrix (mask_fixels); + // Also pre-compute fixel-fixel weights for smoothing. + Stats::CFE::norm_connectivity_matrix_type smoothing_weights (mask_fixels); bool do_smoothing = false; const float gaussian_const2 = 2.0 * smooth_std_dev * smooth_std_dev; @@ -307,7 +308,7 @@ void run() { { ProgressBar progress ("normalising and thresholding fixel-fixel connectivity matrix", num_fixels); - for (uint32_t fixel = 0; fixel < num_fixels; ++fixel) { + for (index_type fixel = 0; fixel < num_fixels; ++fixel) { mask.index(0) = fixel; const int32_t row = fixel2row[fixel]; @@ -316,7 +317,6 @@ void run() { // Here, the connectivity matrix needs to be modified to reflect the // fact that fixel indices in the template fixel image may not // correspond to rows in the statistical analysis - std::map new_connectivity_matrix_row; connectivity_value_type sum_weights = 0.0; for (auto& it : connectivity_matrix[fixel]) { @@ -332,34 +332,32 @@ void run() { const value_type distance = std::sqrt (Math::pow2 (positions[fixel][0] - positions[it.first][0]) + Math::pow2 (positions[fixel][1] - positions[it.first][1]) + Math::pow2 (positions[fixel][2] - positions[it.first][2])); - const connectivity_value_type smoothing_weight = connectivity * gaussian_const1 * std::exp (-std::pow (distance, 2) / gaussian_const2); + const connectivity_value_type smoothing_weight = connectivity * gaussian_const1 * std::exp (-Math::pow2 (distance) / gaussian_const2); if (smoothing_weight >= connectivity_threshold) { - smoothing_weights[row][fixel2row[it.first]] = smoothing_weight; + smoothing_weights[row].push_back (Stats::CFE::NormMatrixElement (fixel2row[it.first], smoothing_weight)); sum_weights += smoothing_weight; } } // Here we pre-exponentiate each connectivity value by C - new_connectivity_matrix_row[fixel2row[it.first]] = std::pow (connectivity, cfe_c); + norm_connectivity_matrix[row].push_back (Stats::CFE::NormMatrixElement (fixel2row[it.first], std::pow (connectivity, cfe_c))); } } // Make sure the fixel is fully connected to itself - new_connectivity_matrix_row[row] = Stats::CFE::connectivity (1.0); - smoothing_weights[row][row] = connectivity_value_type(gaussian_const1); + norm_connectivity_matrix[row].push_back (Stats::CFE::NormMatrixElement (uint32_t(row), connectivity_value_type(1.0))); + smoothing_weights[row].push_back (Stats::CFE::NormMatrixElement (uint32_t(row), connectivity_value_type(gaussian_const1))); sum_weights += connectivity_value_type(gaussian_const1); // Normalise smoothing weights const connectivity_value_type norm_factor = connectivity_value_type(1.0) / sum_weights; for (auto i : smoothing_weights[row]) - i.second *= norm_factor; - - // Write the new connectivity matrix data to the relevant structure - // Note that due to fixel masking, this may not be the same row as - // the fixel index - connectivity_matrix[row] = std::move (new_connectivity_matrix_row); - // Force deallocation of memory used - std::map().swap (new_connectivity_matrix_row); + i.normalise (norm_factor); + + // Force deallocation of memory used for this fixel in the original matrix + std::map().swap (connectivity_matrix[fixel]); + } else { + // If fixel is not in the mask, tract_processor should never assign // any connections to it assert (connectivity_matrix[fixel].empty()); @@ -370,9 +368,11 @@ void run() { } } - // If the connectivity matrix has shrunk in size due to the masking of fixels, - // throw out any now-erroneous data off the end of the vector - connectivity_matrix.resize (mask_fixels); + // The connectivity matrix is now in vector rather than matrix form; + // throw out the structure holding the original data + // (Note however that all entries in the original structure should + // have been deleted during the prior loop) + Stats::CFE::init_connectivity_matrix_type().swap (connectivity_matrix); Header output_header (header); @@ -390,7 +390,7 @@ void run() { matrix_type data (mask_fixels, identifiers.size()); data.setZero(); { - ProgressBar progress ("loading input images", identifiers.size()); + ProgressBar progress (std::string ("loading input images") + (do_smoothing ? " and smoothing" : ""), identifiers.size()); for (size_t subject = 0; subject < identifiers.size(); subject++) { LogLevelLatch log_level (0); @@ -415,7 +415,7 @@ void run() { for (size_t fixel = 0; fixel < mask_fixels; ++fixel) { value_type value = 0.0; for (auto i : smoothing_weights[fixel]) - value += subject_data_vector[i.first] * i.second; + value += subject_data_vector[i.index()] * i.value(); data (fixel, subject) = value; } } else { @@ -425,6 +425,8 @@ void run() { } } + // Free the memory occupied by the data smoothing filter; no longer required + Stats::CFE::norm_connectivity_matrix_type().swap (smoothing_weights); if (!data.allFinite()) throw Exception ("input data contains non-finite value(s)"); @@ -446,7 +448,7 @@ void run() { Math::Stats::GLMTTest glm_ttest (data, design, contrast); std::shared_ptr cfe_integrator; - cfe_integrator.reset (new Stats::CFE::Enhancer (connectivity_matrix, cfe_dh, cfe_e, cfe_h)); + cfe_integrator.reset (new Stats::CFE::Enhancer (norm_connectivity_matrix, cfe_dh, cfe_e, cfe_h)); vector_type empirical_cfe_statistic; // If performing non-stationarity adjustment we need to pre-compute the empirical CFE statistic diff --git a/src/stats/cfe.cpp b/src/stats/cfe.cpp index 1f012fda53..83d21dbef5 100644 --- a/src/stats/cfe.cpp +++ b/src/stats/cfe.cpp @@ -23,11 +23,11 @@ namespace MR - TrackProcessor::TrackProcessor (Image& fixel_indexer, + TrackProcessor::TrackProcessor (Image& fixel_indexer, const vector& fixel_directions, Image& fixel_mask, vector& fixel_TDI, - vector >& connectivity_matrix, + init_connectivity_matrix_type& connectivity_matrix, const value_type angular_threshold) : fixel_indexer (fixel_indexer) , fixel_directions (fixel_directions), @@ -41,21 +41,21 @@ namespace MR bool TrackProcessor::operator() (const SetVoxelDir& in) { // For each voxel tract tangent, assign to a fixel - vector tract_fixel_indices; + vector tract_fixel_indices; for (SetVoxelDir::const_iterator i = in.begin(); i != in.end(); ++i) { assign_pos_of (*i).to (fixel_indexer); fixel_indexer.index(3) = 0; - uint32_t num_fixels = fixel_indexer.value(); + const index_type num_fixels = fixel_indexer.value(); if (num_fixels > 0) { fixel_indexer.index(3) = 1; - uint32_t first_index = fixel_indexer.value(); - uint32_t last_index = first_index + num_fixels; + const index_type first_index = fixel_indexer.value(); + const index_type last_index = first_index + num_fixels; // Note: Streamlines can still be assigned to a fixel that is outside the mask; // however this will not be permitted to contribute to the matrix - uint32_t closest_fixel_index = num_fixels; + index_type closest_fixel_index = num_fixels; value_type largest_dp = 0.0; const direction_type dir (i->get_dir().normalized()); - for (uint32_t j = first_index; j < last_index; ++j) { + for (index_type j = first_index; j < last_index; ++j) { const value_type dp = std::abs (dir.dot (fixel_directions[j])); if (dp > largest_dp) { largest_dp = dp; @@ -92,11 +92,11 @@ namespace MR - Enhancer::Enhancer (const vector >& connectivity_map, + Enhancer::Enhancer (const norm_connectivity_matrix_type& connectivity_matrix, const value_type dh, const value_type E, const value_type H) : - connectivity_map (connectivity_map), + connectivity_matrix (connectivity_matrix), dh (dh), E (E), H (H) { } @@ -107,13 +107,13 @@ namespace MR { enhanced_stats = vector_type::Zero (stats.size()); value_type max_enhanced_stat = 0.0; - for (size_t fixel = 0; fixel < connectivity_map.size(); ++fixel) { - std::map::const_iterator connected_fixel; + vector::const_iterator connected_fixel; + for (size_t fixel = 0; fixel < connectivity_matrix.size(); ++fixel) { for (value_type h = this->dh; h < stats[fixel]; h += this->dh) { value_type extent = 0.0; - for (connected_fixel = connectivity_map[fixel].begin(); connected_fixel != connectivity_map[fixel].end(); ++connected_fixel) - if (stats[connected_fixel->first] > h) - extent += connected_fixel->second.value; + for (connected_fixel = connectivity_matrix[fixel].begin(); connected_fixel != connectivity_matrix[fixel].end(); ++connected_fixel) + if (stats[connected_fixel->index()] > h) + extent += connected_fixel->value(); enhanced_stats[fixel] += std::pow (extent, E) * std::pow (h, H); } if (enhanced_stats[fixel] > max_enhanced_stat) diff --git a/src/stats/cfe.h b/src/stats/cfe.h index 1130be419e..6bf5ebd99e 100644 --- a/src/stats/cfe.h +++ b/src/stats/cfe.h @@ -17,6 +17,7 @@ #include "image.h" #include "image_helpers.h" +#include "types.h" #include "math/math.h" #include "math/stats/typedefs.h" @@ -30,11 +31,11 @@ namespace MR namespace CFE { + using index_type = uint32_t; using value_type = Math::Stats::value_type; using vector_type = Math::Stats::vector_type; using connectivity_value_type = float; using direction_type = Eigen::Matrix; - using connectivity_vector_type = Eigen::Array; using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir; @@ -51,6 +52,30 @@ namespace MR }; + // A class to store fixel index / connectivity value pairs + // only after the connectivity matrix has been thresholded / normalised + class NormMatrixElement + { + public: + NormMatrixElement (const index_type fixel_index, + const connectivity_value_type connectivity_value) : + fixel_index (fixel_index), + connectivity_value (connectivity_value) { } + index_type index() const { return fixel_index; } + connectivity_value_type value() const { return connectivity_value; } + void normalise (const connectivity_value_type norm_factor) { connectivity_value *= norm_factor; } + private: + const index_type fixel_index; + connectivity_value_type connectivity_value; + }; + + + + // Different types are used depending on whether the connectivity matrix + // is in the process of being built, or whether it has been normalised + using init_connectivity_matrix_type = vector>; + using norm_connectivity_matrix_type = vector>; + /** @@ -59,21 +84,21 @@ namespace MR class TrackProcessor { MEMALIGN(TrackProcessor) public: - TrackProcessor (Image& fixel_indexer, + TrackProcessor (Image& fixel_indexer, const vector& fixel_directions, Image& fixel_mask, vector& fixel_TDI, - vector >& connectivity_matrix, + init_connectivity_matrix_type& connectivity_matrix, const value_type angular_threshold); bool operator () (const SetVoxelDir& in); private: - Image fixel_indexer; + Image fixel_indexer; const vector& fixel_directions; Image fixel_mask; vector& fixel_TDI; - vector >& connectivity_matrix; + init_connectivity_matrix_type& connectivity_matrix; const value_type angular_threshold_dp; }; @@ -82,7 +107,7 @@ namespace MR class Enhancer : public Stats::EnhancerBase { MEMALIGN (Enhancer) public: - Enhancer (const vector >& connectivity_map, + Enhancer (const norm_connectivity_matrix_type& connectivity_matrix, const value_type dh, const value_type E, const value_type H); @@ -90,7 +115,7 @@ namespace MR protected: - const vector >& connectivity_map; + const norm_connectivity_matrix_type& connectivity_matrix; const value_type dh, E, H; }; From 3769f48b80898495245feab934db9c004f35bcf1 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 4 Jul 2017 12:15:51 +1000 Subject: [PATCH 06/11] fixelcfestats: Make sure connectivity matrix access functions are inline --- lib/mrtrix3/app.py | 4 +--- src/stats/cfe.h | 6 +++--- 2 files changed, 4 insertions(+), 6 deletions(-) diff --git a/lib/mrtrix3/app.py b/lib/mrtrix3/app.py index ae1306e8c8..d79964beb6 100644 --- a/lib/mrtrix3/app.py +++ b/lib/mrtrix3/app.py @@ -7,13 +7,11 @@ # - cmdline.addCitation(), cmdline.addDescription(), cmdline.setCopyright() as needed # - Add arguments and options to 'cmdline' as needed # - parse() -# - checkOutputFile() as needed +# - checkOutputPath() as needed # - makeTempDir() if the script requires a temporary directory # - gotoTempDir() if the script is using a temporary directory # - complete() # -# checkOutputFile() can be called at any time after parse() to check if an intended output -# file already exists diff --git a/src/stats/cfe.h b/src/stats/cfe.h index 6bf5ebd99e..f6029d1a89 100644 --- a/src/stats/cfe.h +++ b/src/stats/cfe.h @@ -61,9 +61,9 @@ namespace MR const connectivity_value_type connectivity_value) : fixel_index (fixel_index), connectivity_value (connectivity_value) { } - index_type index() const { return fixel_index; } - connectivity_value_type value() const { return connectivity_value; } - void normalise (const connectivity_value_type norm_factor) { connectivity_value *= norm_factor; } + FORCE_INLINE index_type index() const { return fixel_index; } + FORCE_INLINE connectivity_value_type value() const { return connectivity_value; } + FORCE_INLINE void normalise (const connectivity_value_type norm_factor) { connectivity_value *= norm_factor; } private: const index_type fixel_index; connectivity_value_type connectivity_value; From 840beacf53872a2e9e1eff53211ee6f8eba3e810 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 4 Jul 2017 12:55:00 +1000 Subject: [PATCH 07/11] fixelcfestats: Minor tweaks --- cmd/fixelcfestats.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index ac06fb7ac6..14adb18c98 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -173,7 +173,7 @@ void run() { mask_fixels = 0; for (mask.index(0) = 0; mask.index(0) != num_fixels; ++mask.index(0)) fixel2row[mask.index(0)] = mask.value() ? mask_fixels++ : -1; - CONSOLE ("Mask contains " + str(mask_fixels) + " fixels"); + CONSOLE ("Fixel mask contains " + str(mask_fixels) + " fixels"); } else { Header data_header; data_header.ndim() = 3; @@ -183,10 +183,12 @@ void run() { data_header.spacing(0) = data_header.spacing(1) = data_header.spacing(2) = NaN; data_header.stride(0) = 1; data_header.stride(1) = 2; data_header.stride(2) = 3; mask = Image::scratch (data_header, "scratch fixel mask"); - for (mask.index(0) = 0; mask.index(0) != num_fixels; ++mask.index(0)) + for (index_type f = 0; f != num_fixels; ++f) { + mask.index(0) = f; mask.value() = true; + } mask_fixels = num_fixels; - for (uint32_t f = 0; f != num_fixels; ++f) + for (index_type f = 0; f != num_fixels; ++f) fixel2row[f] = f; } @@ -306,6 +308,14 @@ void run() { gaussian_const1 = 1.0 / (smooth_std_dev * std::sqrt (2.0 * Math::pi)); } + for (size_t f = 0; f != num_fixels; ++f) { + mask.index(0) = f; + if (!mask.value()) { + TRACE; + throw Exception ("False value in mask"); + } + } + { ProgressBar progress ("normalising and thresholding fixel-fixel connectivity matrix", num_fixels); for (index_type fixel = 0; fixel < num_fixels; ++fixel) { @@ -493,6 +503,9 @@ void run() { uncorrected_pvalues_neg.reset (new vector_type (mask_fixels)); } + // FIXME fixelcfestats is hanging here for some reason... + // Even when no mask is supplied + if (permutations.size()) { Stats::PermTest::run_permutations (permutations, glm_ttest, cfe_integrator, empirical_cfe_statistic, cfe_output, cfe_output_neg, From f7460550e0dfbed46775f0131c3ce88916a3647c Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 5 Jul 2017 13:49:12 +1000 Subject: [PATCH 08/11] mrview Fixel plot: Compatibility with fixel data files containing NaNs --- src/gui/mrview/tool/fixel/base_fixel.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/gui/mrview/tool/fixel/base_fixel.cpp b/src/gui/mrview/tool/fixel/base_fixel.cpp index 0084937d0d..a21333dae5 100644 --- a/src/gui/mrview/tool/fixel/base_fixel.cpp +++ b/src/gui/mrview/tool/fixel/base_fixel.cpp @@ -118,35 +118,37 @@ namespace MR "void main() {\n"; if (fixel.use_discard_lower()) - source += " if (v_threshold[0] < lower) return;\n"; + source += " if (v_threshold[0] < lower || isnan(v_threshold[0])) return;\n"; if (fixel.use_discard_upper()) - source += " if (v_threshold[0] > upper) return;\n"; + source += " if (v_threshold[0] > upper || isnan(v_threshold[0])) return;\n"; switch (scale_type) { case Unity: - source += " vec4 line_offset = length_mult * vec4 (v_dir[0], 0);\n"; + source += " vec4 line_offset = length_mult * vec4 (v_dir[0], 0);\n"; break; case Value: - source += " vec4 line_offset = length_mult * v_scale[0] * vec4 (v_dir[0], 0);\n"; + source += " if (isnan(v_scale[0])) return;\n" + " vec4 line_offset = length_mult * v_scale[0] * vec4 (v_dir[0], 0);\n"; break; } switch (color_type) { case CValue: if (!ColourMap::maps[colourmap].special) { - source += " float amplitude = clamp ("; + source += " if (isnan(v_colour[0])) return;\n" + " float amplitude = clamp ("; if (fixel.scale_inverted()) source += "1.0 -"; source += " scale * (v_colour[0] - offset), 0.0, 1.0);\n"; } source += - std::string (" vec3 color;\n") + + std::string (" vec3 color;\n") + ColourMap::maps[colourmap].glsl_mapping + - " fColour = color;\n"; + " fColour = color;\n"; break; case Direction: source += - " fColour = normalize (abs (v_dir[0]));\n"; + " fColour = normalize (abs (v_dir[0]));\n"; break; default: break; From 07289db7c101a4d4f2bb9e2a66cf0163629de218 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 5 Jul 2017 13:54:04 +1000 Subject: [PATCH 09/11] Stats::CFE::NormMatrixElement class: Eigen 3.3 compatibility --- src/stats/cfe.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stats/cfe.h b/src/stats/cfe.h index f6029d1a89..5c62d8cb7d 100644 --- a/src/stats/cfe.h +++ b/src/stats/cfe.h @@ -55,7 +55,7 @@ namespace MR // A class to store fixel index / connectivity value pairs // only after the connectivity matrix has been thresholded / normalised class NormMatrixElement - { + { NOMEMALIGN public: NormMatrixElement (const index_type fixel_index, const connectivity_value_type connectivity_value) : From 0f13690f74a0abab0a995151118ca52ba5d97571 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Wed, 5 Jul 2017 16:33:51 +1000 Subject: [PATCH 10/11] fixelcfestats: Remove vestigial debugging call --- cmd/fixelcfestats.cpp | 8 -------- docs/reference/commands/fixelcfestats.rst | 2 +- 2 files changed, 1 insertion(+), 9 deletions(-) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 14adb18c98..800c8a6bdd 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -308,14 +308,6 @@ void run() { gaussian_const1 = 1.0 / (smooth_std_dev * std::sqrt (2.0 * Math::pi)); } - for (size_t f = 0; f != num_fixels; ++f) { - mask.index(0) = f; - if (!mask.value()) { - TRACE; - throw Exception ("False value in mask"); - } - } - { ProgressBar progress ("normalising and thresholding fixel-fixel connectivity matrix", num_fixels); for (index_type fixel = 0; fixel < num_fixels; ++fixel) { diff --git a/docs/reference/commands/fixelcfestats.rst b/docs/reference/commands/fixelcfestats.rst index 209c92d544..da3d415250 100644 --- a/docs/reference/commands/fixelcfestats.rst +++ b/docs/reference/commands/fixelcfestats.rst @@ -94,7 +94,7 @@ Raffelt, D.; Smith, RE.; Ridgway, GR.; Tournier, JD.; Vaughan, DN.; Rose, S.; He -**Author:** David Raffelt (david.raffelt@florey.edu.au) +**Author:** David Raffelt (david.raffelt@florey.edu.au) and Robert E. Smith (robert.smith@florey.edu.au) **Copyright:** Copyright (c) 2008-2017 the MRtrix3 contributors. From 9866913d476a415d6a4d30b3c9cdceb6deabdfd0 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Fri, 7 Jul 2017 17:45:48 +1000 Subject: [PATCH 11/11] fixelcfestats: Add doc information regarding output when using -mask option --- cmd/fixelcfestats.cpp | 7 +++++++ docs/reference/commands/fixelcfestats.rst | 5 +++++ 2 files changed, 12 insertions(+) diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 800c8a6bdd..fbcd485a22 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -55,6 +55,13 @@ void usage () SYNOPSIS = "Fixel-based analysis using connectivity-based fixel enhancement and non-parametric permutation testing"; + DESCRIPTION + + "Note that if the -mask option is used, the output fixel directory will still contain the same set of fixels as that " + "present in the input fixel template, in order to retain fixel correspondence. However a consequence of this is that " + "all fixels in the template will be initialy visible when the output fixel directory is loaded in mrview. Those fixels " + "outside the processing mask will immediately disappear from view as soon as any data-file-based fixel colouring or " + "thresholding is applied."; + REFERENCES + "Raffelt, D.; Smith, RE.; Ridgway, GR.; Tournier, JD.; Vaughan, DN.; Rose, S.; Henderson, R.; Connelly, A." // Internal "Connectivity-based fixel enhancement: Whole-brain statistical analysis of diffusion MRI measures in the presence of crossing fibres. \n" diff --git a/docs/reference/commands/fixelcfestats.rst b/docs/reference/commands/fixelcfestats.rst index da3d415250..7e8655a4bd 100644 --- a/docs/reference/commands/fixelcfestats.rst +++ b/docs/reference/commands/fixelcfestats.rst @@ -22,6 +22,11 @@ Usage - *tracks*: the tracks used to determine fixel-fixel connectivity - *out_fixel_directory*: the output directory where results will be saved. Will be created if it does not exist +Description +----------- + +Note that if the -mask option is used, the output fixel directory will still contain the same set of fixels as that present in the input fixel template, in order to retain fixel correspondence. However a consequence of this is that all fixels in the template will be initialy visible when the output fixel directory is loaded in mrview. Those fixels outside the processing mask will immediately disappear from view as soon as any data-file-based fixel colouring or thresholding is applied. + Options -------