Skip to content

Commit

Permalink
fixel2tsf: Don't propagate non-finite values
Browse files Browse the repository at this point in the history 
With the introduction of the -mask option in fixelcfestats (#1030), fixel data files may contain the value NaN to indicate fixels that were outside the processing mask. However if such values were read within fixel2tsf, and hence propagated to the output track scalar file, these would be erroneously interpreted as track delimiters. This fix replaces non-finite values with 0.0 in fixel2tsf to overcome this limitation of the .tck / .tsf formats.
Some assert() calls have also been added to the relevant code in mrview to ensure that any loaded .tsf file does indeed align perfectly with the corresponding track file.
  • Loading branch information
@Lestropie
Lestropie committed Dec 14, 2017
1 parent 74c1436 commit b2ff414
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 27 deletions.
14 changes: 9 additions & 5 deletions cmd/fixel2tsf.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,15 +98,15 @@ void run ()
while (reader (tck)) {
SetVoxelDir dixels;
mapper (tck, dixels);
vector<float> scalars (tck.size(), 0.0);
vector<float> scalars (tck.size(), 0.0f);
for (size_t p = 0; p < tck.size(); ++p) {
voxel_pos = transform.scanner2voxel * tck[p].cast<default_type> ();
for (SetVoxelDir::const_iterator d = dixels.begin(); d != dixels.end(); ++d) {
if ((int)round(voxel_pos[0]) == (*d)[0] && (int)round(voxel_pos[1]) == (*d)[1] && (int)round(voxel_pos[2]) == (*d)[2]) {
assign_pos_of (*d).to (in_index_image);
Eigen::Vector3f dir = d->get_dir().cast<float>();
dir.normalize();
float largest_dp = 0.0;
float largest_dp = 0.0f;
int32_t closest_fixel_index = -1;

in_index_image.index(3) = 0;
Expand All @@ -116,17 +116,21 @@ void run ()

for (size_t fixel = 0; fixel < num_fixels_in_voxel; ++fixel) {
in_directions_image.index(0) = offset + fixel;
float dp = std::abs (dir.dot (Eigen::Vector3f (in_directions_image.row(1))));
const float dp = std::abs (dir.dot (Eigen::Vector3f (in_directions_image.row(1))));
if (dp > largest_dp) {
largest_dp = dp;
closest_fixel_index = fixel;
}
}
if (largest_dp > angular_threshold_dp) {
in_data_image.index(0) = offset + closest_fixel_index;
scalars[p] = in_data_image.value();
const float value = in_data_image.value();
if (std::isfinite (value))
scalars[p] = in_data_image.value();
else
scalars[p] = 0.0f;
} else {
scalars[p] = 0.0;
scalars[p] = 0.0f;
}
break;
}
Expand Down
56 changes: 37 additions & 19 deletions src/gui/mrview/tool/tractography/tractogram.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -616,8 +616,8 @@ namespace MR

while (file (tck)) {

size_t N = tck.size();
if(!N) continue;
const size_t N = tck.size();
if (!N) continue;

// Pre padding
// To support downsampling, we want to ensure that the starting track vertex
Expand All @@ -633,7 +633,7 @@ namespace MR
// Similarly, to support downsampling, we also want to ensure the final track vertex
// will be used even we're using a stride > 1
for (size_t i = 0; i < track_padding; ++i)
buffer.push_back(tck.back());
buffer.push_back (tck.back());

sizes.push_back (N);
tck_count++;
Expand All @@ -642,9 +642,8 @@ namespace MR

endpoint_tangents.push_back ((tck.back() - tck.front()).normalized());
}
if (buffer.size()) {
if (buffer.size())
load_tracks_onto_GPU (buffer, starts, sizes, tck_count);
}
file.close();
ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;
}
Expand Down Expand Up @@ -682,6 +681,7 @@ namespace MR
}
load_end_colours_onto_GPU (buffer);
}
assert (colour_buffers.size() == vertex_buffers.size());
// Don't need this now that we've initialised the GPU buffers
endpoint_tangents.clear();
ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;
Expand All @@ -708,11 +708,13 @@ namespace MR
DWI::Tractography::Properties scalar_properties;
DWI::Tractography::ScalarReader<float> file (filename, scalar_properties);
DWI::Tractography::check_properties_match (properties, scalar_properties, ".tck / .tsf");
size_t tck_count = 0;
while (file (tck_scalar)) {

size_t tck_size = tck_scalar.size();
const size_t tck_size = tck_scalar.size();
assert (tck_size == size_t(track_sizes[intensity_scalar_buffers.size()][tck_count]));

if(!tck_size)
if (!tck_size)
continue;

// Pre padding to coincide with tracks buffer
Expand All @@ -729,11 +731,13 @@ namespace MR
for (size_t i = 0; i < track_padding; ++i)
buffer.push_back (tck_scalar.back());

++tck_count;

if (buffer.size() >= MAX_BUFFER_SIZE)
load_intensity_scalars_onto_GPU (buffer);
load_intensity_scalars_onto_GPU (buffer, tck_count);
}
if (buffer.size())
load_intensity_scalars_onto_GPU (buffer);
load_intensity_scalars_onto_GPU (buffer, tck_count);
file.close();
} else {
const Eigen::VectorXf scalars = MR::load_vector<float> (filename);
Expand All @@ -745,7 +749,8 @@ namespace MR
size_t running_index = 0;

for (size_t buffer_index = 0; buffer_index != vertex_buffers.size(); ++buffer_index) {
const size_t num_tracks = num_tracks_per_buffer[buffer_index];

size_t num_tracks = num_tracks_per_buffer[buffer_index];
vector<GLint>& track_lengths (original_track_sizes[buffer_index]);

for (size_t index = 0; index != num_tracks; ++index, ++running_index) {
Expand All @@ -766,9 +771,10 @@ namespace MR
value_min = std::min (value_min, value);
}

load_intensity_scalars_onto_GPU (buffer);
load_intensity_scalars_onto_GPU (buffer, num_tracks);
}
}
assert (intensity_scalar_buffers.size() == vertex_buffers.size());
intensity_scalar_filename = filename;
this->set_windowing (value_min, value_max);
if (!std::isfinite (greaterthan))
Expand Down Expand Up @@ -797,11 +803,13 @@ namespace MR
DWI::Tractography::Properties scalar_properties;
DWI::Tractography::ScalarReader<float> file (filename, scalar_properties);
DWI::Tractography::check_properties_match (properties, scalar_properties, ".tck / .tsf");
size_t tck_count = 0;
while (file (tck_scalar)) {

size_t tck_size = tck_scalar.size();
const size_t tck_size = tck_scalar.size();
assert (tck_size == size_t(track_sizes[intensity_scalar_buffers.size()][tck_count]));

if(!tck_size)
if (!tck_size)
continue;

// Pre padding to coincide with tracks buffer
Expand All @@ -818,11 +826,13 @@ namespace MR
for (size_t i = 0; i < track_padding; ++i)
buffer.push_back (tck_scalar.back());

++tck_count;

if (buffer.size() >= MAX_BUFFER_SIZE)
load_threshold_scalars_onto_GPU (buffer);
load_threshold_scalars_onto_GPU (buffer, tck_count);
}
if (buffer.size())
load_threshold_scalars_onto_GPU (buffer);
load_threshold_scalars_onto_GPU (buffer, tck_count);
file.close();
} else {
const Eigen::VectorXf scalars = MR::load_vector<float> (filename);
Expand All @@ -834,7 +844,8 @@ namespace MR
size_t running_index = 0;

for (size_t buffer_index = 0; buffer_index != vertex_buffers.size(); ++buffer_index) {
const size_t num_tracks = num_tracks_per_buffer[buffer_index];

size_t num_tracks = num_tracks_per_buffer[buffer_index];
vector<GLint>& track_lengths (original_track_sizes[buffer_index]);

for (size_t index = 0; index != num_tracks; ++index, ++running_index) {
Expand All @@ -855,9 +866,10 @@ namespace MR
threshold_min = std::min (threshold_min, value);
}

load_threshold_scalars_onto_GPU (buffer);
load_threshold_scalars_onto_GPU (buffer, num_tracks);
}
}
assert (threshold_scalar_buffers.size() == vertex_buffers.size());
threshold_scalar_filename = filename;
greaterthan = threshold_max;
lessthan = threshold_min;
Expand Down Expand Up @@ -996,10 +1008,12 @@ namespace MR



void Tractogram::load_intensity_scalars_onto_GPU (vector<float>& buffer)
void Tractogram::load_intensity_scalars_onto_GPU (vector<float>& buffer, size_t& tck_count)
{
ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;

assert (num_tracks_per_buffer[intensity_scalar_buffers.size()] == tck_count);

GLuint vertexbuffer;
gl::GenBuffers (1, &vertexbuffer);
gl::BindBuffer (gl::ARRAY_BUFFER, vertexbuffer);
Expand All @@ -1009,17 +1023,20 @@ namespace MR

intensity_scalar_buffers.push_back (vertexbuffer);
buffer.clear();
tck_count = 0;

ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;
}




void Tractogram::load_threshold_scalars_onto_GPU (vector<float>& buffer)
void Tractogram::load_threshold_scalars_onto_GPU (vector<float>& buffer, size_t& tck_count)
{
ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;

assert (num_tracks_per_buffer[threshold_scalar_buffers.size()] == tck_count);

GLuint vertexbuffer;
gl::GenBuffers (1, &vertexbuffer);
gl::BindBuffer (gl::ARRAY_BUFFER, vertexbuffer);
Expand All @@ -1029,6 +1046,7 @@ namespace MR

threshold_scalar_buffers.push_back (vertexbuffer);
buffer.clear();
tck_count = 0;

ASSERT_GL_MRVIEW_CONTEXT_IS_CURRENT;
}
Expand Down
6 changes: 3 additions & 3 deletions src/gui/mrview/tool/tractography/tractogram.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -159,11 +159,11 @@ namespace MR
vector<GLint>& starts,
vector<GLint>& sizes,
size_t& tck_count);

void load_end_colours_onto_GPU (vector<Eigen::Vector3f>&);

void load_intensity_scalars_onto_GPU (vector<float>& buffer);
void load_threshold_scalars_onto_GPU (vector<float>& buffer);
void load_intensity_scalars_onto_GPU (vector<float>& buffer, size_t& tck_count);
void load_threshold_scalars_onto_GPU (vector<float>& buffer, size_t& tck_count);

void render_streamlines ();

Expand Down

0 comments on commit b2ff414

Please sign in to comment.