diff --git a/Modules/DiffusionCmdApps/FiberQuantification/CMakeLists.txt b/Modules/DiffusionCmdApps/FiberQuantification/CMakeLists.txt
index 3782c75..aa0841d 100644
--- a/Modules/DiffusionCmdApps/FiberQuantification/CMakeLists.txt
+++ b/Modules/DiffusionCmdApps/FiberQuantification/CMakeLists.txt
@@ -1,40 +1,41 @@
option(BUILD_DiffusionFiberQuantificationCmdApps "Build commandline tools for diffusion fiber processing" OFF)
if(BUILD_DiffusionFiberQuantificationCmdApps OR MITK_BUILD_ALL_APPS)
# needed include directories
include_directories(
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_CURRENT_BINARY_DIR}
)
# list of diffusion cmdapps
# if an app requires additional dependencies
# they are added after a "^^" and separated by "_"
set( diffusionFiberProcessingcmdapps
TractDensity^^
PrintFiberStatistics^^
DistanceToSegmentation^^
+ Tractometry^^
)
foreach(diffusionFiberProcessingcmdapp ${diffusionFiberProcessingcmdapps})
# extract cmd app name and dependencies
string(REPLACE "^^" "\\;" cmdapp_info ${diffusionFiberProcessingcmdapp})
set(cmdapp_info_list ${cmdapp_info})
list(GET cmdapp_info_list 0 appname)
list(GET cmdapp_info_list 1 raw_dependencies)
string(REPLACE "_" "\\;" dependencies "${raw_dependencies}")
set(dependencies_list ${dependencies})
mitkFunctionCreateCommandLineApp(
NAME ${appname}
DEPENDS MitkDiffusionCmdApps MitkMriSimulation MitkFiberTracking ${dependencies_list}
PACKAGE_DEPENDS
)
endforeach()
if(EXECUTABLE_IS_ENABLED)
MITK_INSTALL_TARGETS(EXECUTABLES ${EXECUTABLE_TARGET})
endif()
endif()
diff --git a/Modules/DiffusionCmdApps/FiberQuantification/Tractometry.cpp b/Modules/DiffusionCmdApps/FiberQuantification/Tractometry.cpp
new file mode 100644
index 0000000..632ca43
--- /dev/null
+++ b/Modules/DiffusionCmdApps/FiberQuantification/Tractometry.cpp
@@ -0,0 +1,171 @@
+/*===================================================================
+
+The Medical Imaging Interaction Toolkit (MITK)
+
+Copyright (c) German Cancer Research Center.
+
+All rights reserved.
+
+This software is distributed WITHOUT ANY WARRANTY; without
+even the implied warranty of MERCHANTABILITY or FITNESS FOR
+A PARTICULAR PURPOSE.
+
+See LICENSE.txt or http://www.mitk.org for details.
+
+===================================================================*/
+
+#include <vector>
+#include <iostream>
+#include <fstream>
+#include <algorithm>
+#include <string>
+
+#include <itkImageFileWriter.h>
+#include <itkMetaDataObject.h>
+#include <itkVectorImage.h>
+#include <mitkImageCast.h>
+
+#include <mitkBaseData.h>
+#include <mitkFiberBundle.h>
+#include "mitkDiffusionCommandLineParser.h"
+#include <mitkLexicalCast.h>
+#include <mitkCoreObjectFactory.h>
+#include <mitkIOUtil.h>
+#include <itkTractDensityImageFilter.h>
+#include <itkTractsToFiberEndingsImageFilter.h>
+#include <mitkTractometry.h>
+
+
+mitk::FiberBundle::Pointer LoadFib(std::string filename)
+{
+ std::vector<mitk::BaseData::Pointer> fibInfile = mitk::IOUtil::Load(filename);
+ if( fibInfile.empty() )
+ std::cout << "File " << filename << " could not be read!";
+ mitk::BaseData::Pointer baseData = fibInfile.at(0);
+ return dynamic_cast<mitk::FiberBundle*>(baseData.GetPointer());
+}
+
+/*!
+\brief Modify input tractogram: fiber resampling, compression, pruning and transformation.
+*/
+int main(int argc, char* argv[])
+{
+ mitkDiffusionCommandLineParser parser;
+
+ parser.setTitle("Tract Density");
+ parser.setCategory("Fiber Quantification Methods");
+ parser.setDescription("Generate tract density image, fiber envelope or fiber endpoints image.");
+ parser.setContributor("MIC");
+
+ parser.setArgumentPrefix("--", "-");
+ parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input fiber bundle", us::Any(), false);
+ parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output tracts", us::Any(), false);
+ parser.addArgument("flip", "", mitkDiffusionCommandLineParser::Bool, "flip:", "flip parcellation", us::Any());
+ parser.addArgument("reference_image", "", mitkDiffusionCommandLineParser::String, "Reference image:", "output image will have geometry of this reference image", us::Any(), false);
+ parser.addArgument("num_points", "", mitkDiffusionCommandLineParser::Int, "num_points:", "num_points", 20);
+ parser.addArgument("type", "", mitkDiffusionCommandLineParser::String, "", "CENTROID, STATIC", us::Any(), false);
+
+
+ std::map<std::string, us::Any> parsedArgs = parser.parseArguments(argc, argv);
+ if (parsedArgs.size()==0)
+ return EXIT_FAILURE;
+
+ int num_points = 20;
+ if (parsedArgs.count("num_points"))
+ num_points = us::any_cast<int>(parsedArgs["num_points"]);
+
+ bool flip = false;
+ if (parsedArgs.count("flip"))
+ flip = us::any_cast<bool>(parsedArgs["flip"]);
+
+ MITK_INFO << "Upsampling: " << num_points;
+
+ std::string reference_image = "";
+ if (parsedArgs.count("reference_image"))
+ reference_image = us::any_cast<std::string>(parsedArgs["reference_image"]);
+
+ std::string type = "CENTROID";
+ if (parsedArgs.count("type"))
+ type = us::any_cast<std::string>(parsedArgs["type"]);
+
+ std::string inFileName = us::any_cast<std::string>(parsedArgs["i"]);
+ std::string outFileName = us::any_cast<std::string>(parsedArgs["o"]);
+
+ try
+ {
+ mitk::FiberBundle::Pointer fib = LoadFib(inFileName);
+
+ mitk::Image::Pointer ref_img;
+ ref_img = mitk::IOUtil::Load<mitk::Image>(reference_image);
+
+ mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy();
+ working_fib->ResampleSpline(1.0);
+
+ itk::Image<float, 3>::Pointer itkImage = itk::Image<float, 3>::New();
+ CastToItkImage(ref_img, itkImage);
+
+ if(type == "CENTROID")
+ mitk::Tractometry::NearestCentroidPointTractometry(itkImage, working_fib, num_points, 1, 99999, nullptr, flip);
+ else
+ mitk::Tractometry::StaticResamplingTractometry(itkImage, working_fib, num_points, nullptr, flip);
+
+// std::vector< double > std_values1;
+// std::vector< double > std_values2;
+// std::vector< double > mean_values;
+
+// for (auto row : output)
+// {
+// float mean = row.mean();
+// double stdev = 0;
+
+// for (unsigned int j=0; j<row.size(); ++j)
+// {
+// double diff = mean - row.get(j);
+// diff *= diff;
+// stdev += diff;
+// }
+// stdev /= row.size();
+// stdev = std::sqrt(stdev);
+
+// clipboard_string += boost::lexical_cast<std::string>(mean);
+// clipboard_string += " ";
+// clipboard_string += boost::lexical_cast<std::string>(stdev);
+// clipboard_string += "\n";
+
+// mean_values.push_back(mean);
+// std_values1.push_back(mean + stdev);
+// std_values2.push_back(mean - stdev);
+// }
+// clipboard_string += "\n";
+
+// data.push_back(mean_values);
+// data.push_back(std_values1);
+// data.push_back(std_values2);
+
+// mitk::DataNode::Pointer new_node = mitk::DataNode::New();
+// new_node->SetData(working_fib);
+// new_node->SetName("binned_centroid");
+// new_node->SetVisibility(true);
+// node->SetVisibility(false);
+// GetDataStorage()->Add(new_node, node);
+
+ mitk::IOUtil::Save(working_fib, outFileName );
+
+ }
+ catch (const itk::ExceptionObject& e)
+ {
+ std::cout << e.what();
+ return EXIT_FAILURE;
+ }
+ catch (std::exception& e)
+ {
+ std::cout << e.what();
+ return EXIT_FAILURE;
+ }
+ catch (...)
+ {
+ std::cout << "ERROR!?!";
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
diff --git a/Modules/FiberTracking/Algorithms/mitkTractometry.cpp b/Modules/FiberTracking/Algorithms/mitkTractometry.cpp
index 6f27c46..67f6cd1 100644
--- a/Modules/FiberTracking/Algorithms/mitkTractometry.cpp
+++ b/Modules/FiberTracking/Algorithms/mitkTractometry.cpp
@@ -1,289 +1,330 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include "mitkTractometry.h"
#define _USE_MATH_DEFINES
#include <math.h>
#include <boost/progress.hpp>
#include <vnl/vnl_matrix.h>
#include <mitkTractClusteringFilter.h>
#include <mitkClusteringMetricEuclideanStd.h>
#include <itkTractDensityImageFilter.h>
namespace mitk{
bool Tractometry::Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly)
{
double d_direct = 0;
double d_flipped = 0;
vtkCell* cell1 = polydata1->GetCell(0);
if (ref_poly!=nullptr)
cell1 = ref_poly->GetCell(0);
auto numPoints1 = cell1->GetNumberOfPoints();
vtkPoints* points1 = cell1->GetPoints();
std::vector<itk::Point<double, 3>> ref_points;
for (int j=0; j<numPoints1; ++j)
{
double* p1 = points1->GetPoint(j);
itk::Point<double, 3> itk_p;
itk_p[0] = p1[0];
itk_p[1] = p1[1];
itk_p[2] = p1[2];
ref_points.push_back(itk_p);
}
vtkCell* cell2 = polydata1->GetCell(i);
vtkPoints* points2 = cell2->GetPoints();
for (int j=0; j<numPoints1; ++j)
{
auto p1 = ref_points.at(j);
double* p2 = points2->GetPoint(j);
d_direct += (p1[0]-p2[0])*(p1[0]-p2[0]) + (p1[1]-p2[1])*(p1[1]-p2[1]) + (p1[2]-p2[2])*(p1[2]-p2[2]);
double* p3 = points2->GetPoint(numPoints1-j-1);
d_flipped += (p1[0]-p3[0])*(p1[0]-p3[0]) + (p1[1]-p3[1])*(p1[1]-p3[1]) + (p1[2]-p3[2])*(p1[2]-p3[2]);
}
if (d_direct>d_flipped)
return true;
return false;
}
void Tractometry::ResampleIfNecessary(mitk::FiberBundle::Pointer fib, unsigned int num_points)
{
auto poly = fib->GetFiberPolyData();
bool resample = false;
for (int i=0; i<poly->GetNumberOfCells(); i++)
{
vtkCell* cell = poly->GetCell(i);
if (cell->GetNumberOfPoints()!=num_points)
{
resample = true;
MITK_INFO << "Resampling required!";
break;
}
}
if (resample)
fib->ResampleToNumPoints(num_points);
}
unsigned int Tractometry::EstimateNumSamplingPoints(itk::Image<unsigned char, 3>::Pointer ref_image, mitk::FiberBundle::Pointer fib, unsigned int voxels)
{
auto spacing = ref_image->GetSpacing();
float f = (spacing[0] + spacing[1] + spacing[2])/3;
float num_voxels_passed = 0;
for (unsigned int i=0; i<fib->GetNumFibers(); ++i)
num_voxels_passed += fib->GetFiberLength(i)/f;
num_voxels_passed /= fib->GetNumFibers();
unsigned int parcels = std::ceil(num_voxels_passed/voxels);
MITK_INFO << "Estimated number of sampling points " << parcels;
return parcels;
}
-std::vector<vnl_vector<float>> Tractometry::NearestCentroidPointTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer working_fib, unsigned int num_points, unsigned int max_centroids, float cluster_size, mitk::FiberBundle::Pointer ref_fib)
+std::vector<vnl_vector<float>> Tractometry::NearestCentroidPointTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer working_fib, unsigned int num_points, unsigned int max_centroids, float cluster_size,
+ mitk::FiberBundle::Pointer ref_fib, bool flip_parcellation)
{
vtkSmartPointer< vtkPolyData > polydata = working_fib->GetFiberPolyData();
vtkSmartPointer< vtkPolyData > ref_polydata = nullptr;
if (ref_fib!=nullptr)
{
ResampleIfNecessary(ref_fib, num_points);
ref_polydata = ref_fib->GetFiberPolyData();
}
auto interpolator = itk::LinearInterpolateImageFunction< itk::Image<float, 3>, float >::New();
interpolator->SetInputImage(itkImage);
+ {
+ auto p1 = polydata->GetCell(polydata->GetNumberOfCells()/2)->GetPoints()->GetPoint(0);
+ itk::Point<double, 3> itk_p;
+ itk_p[0] = p1[0];
+ itk_p[1] = p1[1];
+ itk_p[2] = p1[2];
+
+ float pixelValue = mitk::imv::GetImageValue<float, double>(itk_p, false, interpolator);
+ if (pixelValue > 1.5)
+ {
+ flip_parcellation = true;
+ MITK_INFO << "FLIP";
+ }
+ else
+ flip_parcellation = false;
+ }
+
mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New();
lookupTable->SetType(mitk::LookupTable::MULTILABEL);
std::vector<std::vector<float>> output_temp;
for(unsigned int i=0; i<num_points; ++i)
output_temp.push_back({});
// clustering
std::vector< mitk::ClusteringMetric* > metrics;
metrics.push_back({new mitk::ClusteringMetricEuclideanStd()});
mitk::FiberBundle::Pointer resampled = working_fib->GetDeepCopy();
ResampleIfNecessary(resampled, num_points);
std::vector<mitk::FiberBundle::Pointer> centroids;
std::shared_ptr< mitk::TractClusteringFilter > clusterer = std::make_shared<mitk::TractClusteringFilter>();
int c=0;
while (c<30 && (centroids.empty() || centroids.size()>max_centroids))
{
float cs = cluster_size + cluster_size*c*0.2;
float max_d = 0;
int i=1;
std::vector< float > distances;
while (max_d < resampled->GetGeometry()->GetDiagonalLength()/2)
{
distances.push_back(cs*i);
max_d = cs*i;
++i;
}
clusterer->SetDistances(distances);
clusterer->SetTractogram(resampled);
clusterer->SetMetrics(metrics);
clusterer->SetMergeDuplicateThreshold(cs);
clusterer->SetDoResampling(false);
clusterer->SetNumPoints(num_points);
if (c==29)
clusterer->SetMaxClusters(max_centroids);
clusterer->SetMinClusterSize(1);
clusterer->Update();
centroids = clusterer->GetOutCentroids();
++c;
}
for (unsigned int i=0; i<working_fib->GetNumFibers(); ++i)
{
vtkCell* cell = polydata->GetCell(i);
auto numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
vnl_vector<float> values; values.set_size(numPoints);
for (int j=0; j<numPoints; j++)
{
double* p = points->GetPoint(j);
int min_bin = 0;
float d=999999;
for (auto centroid : centroids)
{
auto centroid_polydata = centroid->GetFiberPolyData();
vtkCell* centroid_cell = centroid_polydata->GetCell(0);
auto centroid_numPoints = centroid_cell->GetNumberOfPoints();
vtkPoints* centroid_points = centroid_cell->GetPoints();
bool centroid_flip = Flip(centroid_polydata, 0, ref_polydata);
+ if (flip_parcellation)
+ centroid_flip = !centroid_flip;
for (int bin=0; bin<centroid_numPoints; ++bin)
{
double* centroid_p;
centroid_p = centroid_points->GetPoint(bin);
float temp_d = std::sqrt((p[0]-centroid_p[0])*(p[0]-centroid_p[0]) + (p[1]-centroid_p[1])*(p[1]-centroid_p[1]) + (p[2]-centroid_p[2])*(p[2]-centroid_p[2]));
if (temp_d<d)
{
d = temp_d;
if (centroid_flip)
min_bin = centroid_numPoints-bin-1;
else
min_bin = bin;
}
}
}
double rgb[3] = {0,0,0};
lookupTable->GetColor(min_bin+1, rgb);
working_fib->ColorSinglePoint(i, j, rgb);
double pixelValue = mitk::imv::GetImageValue<float, double>(mitk::imv::GetItkPoint(p), true, interpolator);
output_temp.at(min_bin).push_back(pixelValue);
}
}
std::vector<vnl_vector<float>> output;
for (auto row_v : output_temp)
{
vnl_vector<float> row; row.set_size(row_v.size());
int i = 0;
for (auto v : row_v)
{
row.put(i, v);
++i;
}
output.push_back(row);
}
return output;
}
-vnl_matrix<float> Tractometry::StaticResamplingTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer working_fib, unsigned int num_points, mitk::FiberBundle::Pointer ref_fib)
+vnl_matrix<float> Tractometry::StaticResamplingTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer working_fib, unsigned int num_points, mitk::FiberBundle::Pointer ref_fib, bool flip_parcellation)
{
ResampleIfNecessary(working_fib, num_points);
vtkSmartPointer< vtkPolyData > polydata = working_fib->GetFiberPolyData();
vtkSmartPointer< vtkPolyData > ref_polydata = nullptr;
if (ref_fib!=nullptr)
{
ResampleIfNecessary(ref_fib, num_points);
ref_polydata = ref_fib->GetFiberPolyData();
}
auto interpolator = itk::LinearInterpolateImageFunction< itk::Image<float, 3>, float >::New();
interpolator->SetInputImage(itkImage);
mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New();
lookupTable->SetType(mitk::LookupTable::MULTILABEL);
vnl_matrix<float> output; output.set_size(num_points, working_fib->GetNumFibers());
output.fill(0.0);
+
+ {
+ auto p1 = polydata->GetCell(polydata->GetNumberOfCells()/2)->GetPoints()->GetPoint(0);
+ itk::Point<double, 3> itk_p;
+ itk_p[0] = p1[0];
+ itk_p[1] = p1[1];
+ itk_p[2] = p1[2];
+
+ float pixelValue = mitk::imv::GetImageValue<float, double>(itk_p, false, interpolator);
+ MITK_INFO << "TEST <<< " << pixelValue;
+ if (pixelValue > 1.5)
+ {
+ flip_parcellation = true;
+ MITK_INFO << "FLIP";
+ }
+ else
+ flip_parcellation = false;
+ }
+
for (unsigned int i=0; i<working_fib->GetNumFibers(); ++i)
{
vtkCell* cell = polydata->GetCell(i);
auto numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
bool flip = Flip(polydata, i, ref_polydata);
+ if (flip_parcellation)
+ flip = !flip;
for (int j=0; j<numPoints; j++)
{
double rgb[3] = {0,0,0};
lookupTable->GetColor(j+1, rgb);
double* p;
if (flip)
{
auto p_idx = numPoints - j - 1;
p = points->GetPoint(p_idx);
working_fib->ColorSinglePoint(i, p_idx, rgb);
}
else
{
p = points->GetPoint(j);
working_fib->ColorSinglePoint(i, j, rgb);
}
double pixelValue = mitk::imv::GetImageValue<float, double>(mitk::imv::GetItkPoint(p), true, interpolator);
output.put(j, i, pixelValue);
}
}
return output;
}
}
diff --git a/Modules/FiberTracking/Algorithms/mitkTractometry.h b/Modules/FiberTracking/Algorithms/mitkTractometry.h
index 3e6c64b..2c1f091 100644
--- a/Modules/FiberTracking/Algorithms/mitkTractometry.h
+++ b/Modules/FiberTracking/Algorithms/mitkTractometry.h
@@ -1,59 +1,59 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#ifndef TractometryFilter_h
#define TractometryFilter_h
// MITK
#include <MitkFiberTrackingExports.h>
#include <mitkPlanarEllipse.h>
#include <mitkFiberBundle.h>
#include <mitkClusteringMetric.h>
// ITK
#include <itkProcessObject.h>
// VTK
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkPoints.h>
#include <vtkPolyLine.h>
namespace mitk{
/**
* \brief */
class MITKFIBERTRACKING_EXPORT Tractometry
{
public:
- static vnl_matrix<float> StaticResamplingTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer fib, unsigned int num_points, mitk::FiberBundle::Pointer ref_fib);
+ static vnl_matrix<float> StaticResamplingTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer fib, unsigned int num_points, mitk::FiberBundle::Pointer ref_fib, bool flip_parcellation);
- static std::vector<vnl_vector<float>> NearestCentroidPointTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer fib, unsigned int num_points, unsigned int max_centroids, float cluster_size, mitk::FiberBundle::Pointer ref_fib);
+ static std::vector<vnl_vector<float>> NearestCentroidPointTractometry(itk::Image<float, 3>::Pointer itkImage, mitk::FiberBundle::Pointer fib, unsigned int num_points, unsigned int max_centroids, float cluster_size, mitk::FiberBundle::Pointer ref_fib, bool flip_parcellation);
static unsigned int EstimateNumSamplingPoints(itk::Image<unsigned char, 3>::Pointer ref_image, mitk::FiberBundle::Pointer fib, unsigned int voxels);
static void ResampleIfNecessary(mitk::FiberBundle::Pointer fib, unsigned int num_points);
protected:
static bool Flip(vtkSmartPointer< vtkPolyData > polydata1, int i, vtkSmartPointer< vtkPolyData > ref_poly=nullptr);
};
}
#endif
diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
index 70a1423..dc72718 100644
--- a/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
+++ b/Plugins/org.mitk.gui.qt.diffusionimaging.fiberprocessing/src/internal/QmitkTractometryView.cpp
@@ -1,354 +1,354 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
#include <berryISelectionService.h>
#include <berryIWorkbenchWindow.h>
#include "QmitkTractometryView.h"
#include <mitkNodePredicateDataType.h>
#include <mitkNodePredicateDimension.h>
#include <mitkNodePredicateAnd.h>
#include <mitkImageCast.h>
#include <mitkImageAccessByItk.h>
#include <mitkImage.h>
#include <mitkFiberBundle.h>
#include <mitkDiffusionPropertyHelper.h>
#include <mitkITKImageImport.h>
#include <mitkPixelTypeMultiplex.h>
#include <mitkImagePixelReadAccessor.h>
#include <berryIPreferences.h>
#include <berryWorkbenchPlugin.h>
#include <berryQtPreferences.h>
#include <vtkLookupTable.h>
#include <QClipboard>
#include <mitkLexicalCast.h>
#include <QmitkChartWidget.h>
#include <mitkLookupTable.h>
#include <mitkTractClusteringFilter.h>
#include <mitkClusteringMetricEuclideanStd.h>
#include <itkTractDensityImageFilter.h>
#include <mitkLookupTableProperty.h>
#include <mitkLevelWindowProperty.h>
#include <itkMaskedStatisticsImageFilter.h>
#include <mitkDiffusionFunctionCollection.h>
#include <mitkTractometry.h>
const std::string QmitkTractometryView::VIEW_ID = "org.mitk.views.tractometry";
using namespace mitk;
QmitkTractometryView::QmitkTractometryView()
: QmitkAbstractView()
, m_Controls( nullptr )
, m_Visible(false)
{
}
// Destructor
QmitkTractometryView::~QmitkTractometryView()
{
}
void QmitkTractometryView::CreateQtPartControl( QWidget *parent )
{
// build up qt view, unless already done
if ( !m_Controls )
{
// create GUI widgets from the Qt Designer's .ui file
m_Controls = new Ui::QmitkTractometryViewControls;
m_Controls->setupUi( parent );
connect( m_Controls->m_MethodBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) );
connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartTractometry()) );
mitk::TNodePredicateDataType<mitk::Image>::Pointer imageP = mitk::TNodePredicateDataType<mitk::Image>::New();
mitk::NodePredicateDimension::Pointer dimP = mitk::NodePredicateDimension::New(3);
m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage());
m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(imageP, dimP));
m_Controls->m_ChartWidget->SetXAxisLabel("Tract position");
m_Controls->m_ChartWidget->SetYAxisLabel("Image Value");
m_Controls->m_ChartWidget->SetTheme(QmitkChartWidget::ColorTheme::darkstyle);
}
}
void QmitkTractometryView::SetFocus()
{
}
void QmitkTractometryView::UpdateGui()
{
berry::IWorkbenchPart::Pointer nullPart;
OnSelectionChanged(nullPart, QList<mitk::DataNode::Pointer>(m_CurrentSelection));
}
void QmitkTractometryView::StaticResamplingTractometry(mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector<std::vector<double> > &data, std::string& clipboard_string)
{
itk::Image<float, 3>::Pointer itkImage = itk::Image<float, 3>::New();
CastToItkImage(image, itkImage);
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
unsigned int num_points = m_NumSamplingPoints;
mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy();
- vnl_matrix<float> output = mitk::Tractometry::StaticResamplingTractometry(itkImage, working_fib, num_points, m_ReferenceFib);
+ vnl_matrix<float> output = mitk::Tractometry::StaticResamplingTractometry(itkImage, working_fib, num_points, m_ReferenceFib, false);
std::vector< double > std_values1;
std::vector< double > std_values2;
std::vector< double > mean_values;
for (unsigned int i=0; i<output.rows(); ++i)
{
auto row = output.get_row(i);
float mean = row.mean();
double stdev = 0;
for (unsigned int j=0; j<row.size(); ++j)
{
double diff = mean - row.get(j);
diff *= diff;
stdev += diff;
}
stdev /= row.size();
stdev = std::sqrt(stdev);
clipboard_string += boost::lexical_cast<std::string>(mean);
clipboard_string += " ";
clipboard_string += boost::lexical_cast<std::string>(stdev);
clipboard_string += "\n";
mean_values.push_back(mean);
std_values1.push_back(mean + stdev);
std_values2.push_back(mean - stdev);
}
clipboard_string += "\n";
data.push_back(mean_values);
data.push_back(std_values1);
data.push_back(std_values2);
if (m_Controls->m_ShowBinned->isChecked())
{
mitk::DataNode::Pointer new_node = mitk::DataNode::New();
new_node->SetData(working_fib);
new_node->SetName("binned_static");
new_node->SetVisibility(true);
node->SetVisibility(false);
GetDataStorage()->Add(new_node, node);
}
}
void QmitkTractometryView::NearestCentroidPointTractometry(mitk::Image::Pointer image, mitk::DataNode::Pointer node, std::vector< std::vector< double > >& data, std::string& clipboard_string)
{
mitk::FiberBundle::Pointer fib = dynamic_cast<mitk::FiberBundle*>(node->GetData());
unsigned int num_points = m_NumSamplingPoints;
mitk::FiberBundle::Pointer working_fib = fib->GetDeepCopy();
working_fib->ResampleSpline(1.0);
itk::Image<float, 3>::Pointer itkImage = itk::Image<float, 3>::New();
CastToItkImage(image, itkImage);
- auto output = mitk::Tractometry::NearestCentroidPointTractometry(itkImage, working_fib, num_points, m_Controls->m_MaxCentroids->value(), m_Controls->m_ClusterSize->value(), m_ReferenceFib);
+ auto output = mitk::Tractometry::NearestCentroidPointTractometry(itkImage, working_fib, num_points, m_Controls->m_MaxCentroids->value(), m_Controls->m_ClusterSize->value(), m_ReferenceFib, false);
std::vector< double > std_values1;
std::vector< double > std_values2;
std::vector< double > mean_values;
for (auto row : output)
{
float mean = row.mean();
double stdev = 0;
for (unsigned int j=0; j<row.size(); ++j)
{
double diff = mean - row.get(j);
diff *= diff;
stdev += diff;
}
stdev /= row.size();
stdev = std::sqrt(stdev);
clipboard_string += boost::lexical_cast<std::string>(mean);
clipboard_string += " ";
clipboard_string += boost::lexical_cast<std::string>(stdev);
clipboard_string += "\n";
mean_values.push_back(mean);
std_values1.push_back(mean + stdev);
std_values2.push_back(mean - stdev);
}
clipboard_string += "\n";
data.push_back(mean_values);
data.push_back(std_values1);
data.push_back(std_values2);
if (m_Controls->m_ShowBinned->isChecked())
{
mitk::DataNode::Pointer new_node = mitk::DataNode::New();
new_node->SetData(working_fib);
new_node->SetName("binned_centroid");
new_node->SetVisibility(true);
node->SetVisibility(false);
GetDataStorage()->Add(new_node, node);
}
}
void QmitkTractometryView::Activated()
{
}
void QmitkTractometryView::Deactivated()
{
}
void QmitkTractometryView::Visible()
{
m_Visible = true;
QList<mitk::DataNode::Pointer> selection = GetDataManagerSelection();
berry::IWorkbenchPart::Pointer nullPart;
OnSelectionChanged(nullPart, selection);
}
void QmitkTractometryView::Hidden()
{
m_Visible = false;
}
std::string QmitkTractometryView::RGBToHexString(double *rgb)
{
std::ostringstream os;
for (int i = 0; i < 3; ++i)
{
os << std::setw(2) << std::setfill('0') << std::hex
<< static_cast<int>(rgb[i] * 255);
}
return os.str();
}
void QmitkTractometryView::StartTractometry()
{
m_ReferenceFib = dynamic_cast<mitk::FiberBundle*>(m_CurrentSelection.at(0)->GetData())->GetDeepCopy();
mitk::Image::Pointer image = dynamic_cast<mitk::Image*>(m_Controls->m_ImageBox->GetSelectedNode()->GetData());
MITK_INFO << "Resanmpling reference fibers";
if (m_Controls->m_SamplingPointsBox->value()<3)
{
typedef itk::Image<unsigned char, 3> ParcellationImageType;
ParcellationImageType::Pointer itkImage = ParcellationImageType::New();
CastToItkImage(image, itkImage);
m_NumSamplingPoints = mitk::Tractometry::EstimateNumSamplingPoints(itkImage, m_ReferenceFib, 5);
}
else
m_NumSamplingPoints = m_Controls->m_SamplingPointsBox->value();
m_ReferenceFib->ResampleToNumPoints(m_NumSamplingPoints);
double color[3] = {0,0,0};
mitk::LookupTable::Pointer lookupTable = mitk::LookupTable::New();
lookupTable->SetType(mitk::LookupTable::MULTILABEL);
this->m_Controls->m_ChartWidget->Clear();
std::string clipboardString = "";
int c = 1;
for (auto node : m_CurrentSelection)
{
clipboardString += node->GetName() + "\n";
clipboardString += "mean stdev\n";
std::vector< std::vector< double > > data;
switch (m_Controls->m_MethodBox->currentIndex())
{
case 0:
{
StaticResamplingTractometry( image, node, data, clipboardString );
break;
}
case 1:
{
NearestCentroidPointTractometry( image, node, data, clipboardString );
break;
}
default:
{
StaticResamplingTractometry( image, node, data, clipboardString );
}
}
m_Controls->m_ChartWidget->AddData1D(data.at(0), node->GetName() + " Mean", QmitkChartWidget::ChartType::line);
m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " Mean", QmitkChartWidget::LineStyle::solid);
if (m_Controls->m_StDevBox->isChecked())
{
m_Controls->m_ChartWidget->AddData1D(data.at(1), node->GetName() + " +STDEV", QmitkChartWidget::ChartType::line);
m_Controls->m_ChartWidget->AddData1D(data.at(2), node->GetName() + " -STDEV", QmitkChartWidget::ChartType::line);
m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " +STDEV", QmitkChartWidget::LineStyle::dashed);
m_Controls->m_ChartWidget->SetLineStyle(node->GetName() + " -STDEV", QmitkChartWidget::LineStyle::dashed);
}
lookupTable->GetTableValue(c, color);
this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " Mean", RGBToHexString(color));
if (m_Controls->m_StDevBox->isChecked())
{
color[0] *= 0.5;
color[1] *= 0.5;
color[2] *= 0.5;
this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " +STDEV", RGBToHexString(color));
this->m_Controls->m_ChartWidget->SetColor(node->GetName() + " -STDEV", RGBToHexString(color));
}
this->m_Controls->m_ChartWidget->Show(true);
this->m_Controls->m_ChartWidget->SetShowDataPoints(false);
++c;
}
QApplication::clipboard()->setText(clipboardString.c_str(), QClipboard::Clipboard);
}
void QmitkTractometryView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList<mitk::DataNode::Pointer>& nodes)
{
m_Controls->m_StartButton->setEnabled(false);
if (!m_Visible)
return;
if (m_Controls->m_MethodBox->currentIndex()==0)
m_Controls->m_ClusterFrame->setVisible(false);
else
m_Controls->m_ClusterFrame->setVisible(true);
m_CurrentSelection.clear();
if(m_Controls->m_ImageBox->GetSelectedNode().IsNull())
return;
for (auto node: nodes)
if ( dynamic_cast<mitk::FiberBundle*>(node->GetData()) )
m_CurrentSelection.push_back(node);
if (!m_CurrentSelection.empty())
m_Controls->m_StartButton->setEnabled(true);
}