diff --git a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp
index 0dbc893..7058fb7 100644
--- a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp
+++ b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.cpp
@@ -1,287 +1,310 @@
#include <mitkTrackvis.h>
#include <vtkTransformPolyDataFilter.h>
TrackVisFiberReader::TrackVisFiberReader() { m_Filename = ""; m_FilePointer = nullptr; }
TrackVisFiberReader::~TrackVisFiberReader() { if (m_FilePointer) fclose( m_FilePointer ); }
// Create a TrackVis file and store standard metadata. The file is ready to append fibers.
// ---------------------------------------------------------------------------------------
short TrackVisFiberReader::create(std::string filename , mitk::FiberBundle *fib, bool print_header)
{
m_Header = fib->GetTrackVisHeader();
if (print_header)
this->print_header();
// write the header to the file
m_FilePointer = fopen(filename.c_str(),"w+b");
if (m_FilePointer == nullptr)
{
printf("[ERROR] Unable to create file '%s'\n",filename.c_str());
return 0;
}
if (fwrite((char*)&m_Header, 1, 1000, m_FilePointer) != 1000)
MITK_ERROR << "TrackVis::create : Error occurding during writing fiber.";
this->m_Filename = filename;
return 1;
}
// Open an existing TrackVis file and read metadata information.
// The file pointer is positiond at the beginning of fibers data
// -------------------------------------------------------------
short TrackVisFiberReader::open( std::string filename )
{
m_FilePointer = std::fopen(filename.c_str(), "rb");
if (m_FilePointer == nullptr)
{
printf("[ERROR] Unable to open file '%s'\n",filename.c_str());
return 0;
}
this->m_Filename = filename;
return fread((char*)(&m_Header), 1, 1000, m_FilePointer);
}
short TrackVisFiberReader::write(const mitk::FiberBundle *fib)
{
vtkSmartPointer<vtkPolyData> poly = fib->GetFiberPolyData();
{
mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New();
matrix->Identity();
for (int i=0; i<4; ++i)
for (int j=0; j<4; ++j)
{
if (j<3)
matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]/m_Header.voxel_size[j]);
else
matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]);
}
if (m_Header.voxel_order[0]=='R')
{
matrix->SetElement(0,0,-matrix->GetElement(0,0));
matrix->SetElement(0,1,-matrix->GetElement(0,1));
matrix->SetElement(0,2,-matrix->GetElement(0,2));
matrix->SetElement(0,3,-matrix->GetElement(0,3));
}
if (m_Header.voxel_order[1]=='A')
{
matrix->SetElement(1,0,-matrix->GetElement(1,0));
matrix->SetElement(1,1,-matrix->GetElement(1,1));
matrix->SetElement(1,2,-matrix->GetElement(1,2));
matrix->SetElement(1,3,-matrix->GetElement(1,3));
}
if (m_Header.voxel_order[2]=='I')
{
matrix->SetElement(2,0,-matrix->GetElement(2,0));
matrix->SetElement(2,1,-matrix->GetElement(2,1));
matrix->SetElement(2,2,-matrix->GetElement(2,2));
matrix->SetElement(2,3,-matrix->GetElement(2,3));
}
geometry->SetIndexToWorldTransformByVtkMatrix(matrix);
vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
transformFilter->SetInputData(poly);
transformFilter->SetTransform(geometry->GetVtkTransform()->GetInverse());
transformFilter->Update();
poly = transformFilter->GetOutput();
}
for (unsigned int i=0; i<fib->GetNumFibers(); i++)
{
vtkCell* cell = poly->GetCell(i);
int numPoints = cell->GetNumberOfPoints();
vtkPoints* points = cell->GetPoints();
unsigned int numSaved, pos = 0;
//float* tmp = new float[3*maxSteps];
std::vector< float > tmp;
tmp.reserve(3*numPoints);
numSaved = numPoints;
for(unsigned int i=0; i<numSaved ;i++)
{
double* p = points->GetPoint(i);
// TRK coordinates are corner based, so we have to shift the center based vtk coordinates by half a voxel
tmp[pos++] = p[0] + m_Header.voxel_size[0]/2;
tmp[pos++] = p[1] + m_Header.voxel_size[1]/2;
tmp[pos++] = p[2] + m_Header.voxel_size[2]/2;
}
// write the coordinates to the file
if ( fwrite((char*)&numSaved, 1, 4, m_FilePointer) != 4 )
{
printf( "[ERROR] Problems saving the fiber!\n" );
return 1;
}
if ( fwrite((char*)&(tmp.front()), 1, 4*pos, m_FilePointer) != 4*pos )
{
printf( "[ERROR] Problems saving the fiber!\n" );
return 1;
}
}
return 0;
}
void TrackVisFiberReader::print_header()
{
std::cout << "--------------------------------------------------------" << std::endl;
std::cout << "see http://trackvis.org/docs/?subsect=fileformat" << std::endl;
std::cout << "ONLY vox_to_ras AND voxel_order HEADER ENTRIES ARE USED FOR FIBER COORDINATE TRANSFORMATIONS!" << std::endl;
std::cout << "\nid_string (should be \"TRACK\"): " << m_Header.id_string << std::endl;
std::cout << "dim: [" << std::defaultfloat << m_Header.dim[0] << " " << m_Header.dim[1] << " " << m_Header.dim[2] << "]" << std::endl;
std::cout << "voxel_size: [" << m_Header.voxel_size[0] << " " << m_Header.voxel_size[1] << " " << m_Header.voxel_size[2] << "]" << std::endl;
std::cout << "origin: [" << m_Header.origin[0] << " " << m_Header.origin[1] << " " << m_Header.origin[2] << "]" << std::endl;
std::cout << "vox_to_world: " << std::scientific << std::endl;
std::cout << "[[" << m_Header.vox_to_ras[0][0] << ", " << m_Header.vox_to_ras[0][1] << ", " << m_Header.vox_to_ras[0][2] << ", " << m_Header.vox_to_ras[0][3] << "]" << std::endl;
std::cout << " [" << m_Header.vox_to_ras[1][0] << ", " << m_Header.vox_to_ras[1][1] << ", " << m_Header.vox_to_ras[1][2] << ", " << m_Header.vox_to_ras[1][3] << "]" << std::endl;
std::cout << " [" << m_Header.vox_to_ras[2][0] << ", " << m_Header.vox_to_ras[2][1] << ", " << m_Header.vox_to_ras[2][2] << ", " << m_Header.vox_to_ras[2][3] << "]" << std::endl;
std::cout << " [" << m_Header.vox_to_ras[3][0] << ", " << m_Header.vox_to_ras[3][1] << ", " << m_Header.vox_to_ras[3][2] << ", " << m_Header.vox_to_ras[3][3] << "]]" << std::defaultfloat << std::endl;
std::cout << "voxel_order: " << m_Header.voxel_order[0] << m_Header.voxel_order[1] << m_Header.voxel_order[2] << std::endl;
std::cout << "pad1: " << m_Header.pad1[0] << m_Header.pad1[1] << std::endl;
std::cout << "pad2: " << m_Header.pad2[0] << m_Header.pad2[1] << m_Header.pad2[2] << std::endl;
std::cout << "image_orientation_patient: [" << m_Header.image_orientation_patient[0] << " " << m_Header.image_orientation_patient[1] << " " << m_Header.image_orientation_patient[2] << " " << m_Header.image_orientation_patient[3] << " " << m_Header.image_orientation_patient[4] << " " << m_Header.image_orientation_patient[5] << "]" << std::endl;
std::cout << "invert_x: " << static_cast<bool>(m_Header.invert_x) << std::endl;
std::cout << "invert_y: " << static_cast<bool>(m_Header.invert_y) << std::endl;
std::cout << "invert_z: " << static_cast<bool>(m_Header.invert_z) << std::endl;
std::cout << "swap_xy: " << static_cast<bool>(m_Header.swap_xy) << std::endl;
std::cout << "swap_yz: " << static_cast<bool>(m_Header.swap_yz) << std::endl;
std::cout << "swap_zx: " << static_cast<bool>(m_Header.swap_zx) << std::endl;
std::cout << "n_count: " << m_Header.n_count << std::endl;
std::cout << "version: " << m_Header.version << std::endl;
std::cout << "hdr_size: " << m_Header.hdr_size << std::endl;
std::cout << "\nNot printed: n_scalars, scalar_name, n_properties, property_name, reserved" << std::endl;
std::cout << "--------------------------------------------------------" << std::endl;
}
-short TrackVisFiberReader::read( mitk::FiberBundle* fib, bool use_matrix, bool print_header)
+short TrackVisFiberReader::read( mitk::FiberBundle* fib, bool use_matrix, bool flip_x, bool flip_y, bool flip_z, bool print_header)
{
int numPoints;
vtkSmartPointer<vtkPoints> vtkNewPoints = vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkCellArray> vtkNewCells = vtkSmartPointer<vtkCellArray>::New();
if (print_header)
this->print_header();
while (fread((char*)&numPoints, 1, 4, m_FilePointer)==4)
{
if ( numPoints <= 0 )
{
printf( "[ERROR] Trying to read a fiber with %d points!\n", numPoints );
return -1;
}
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
float tmp[3];
for(int i=0; i<numPoints; i++)
{
if (fread((char*)tmp, 1, 12, m_FilePointer) == 0)
MITK_ERROR << "TrackVis::read: Error during read.";
// TRK coordinates are corner based, so we have to shift them back to center based coordinates
tmp[0] -= m_Header.voxel_size[0]/2;
tmp[1] -= m_Header.voxel_size[1]/2;
tmp[2] -= m_Header.voxel_size[2]/2;
vtkIdType id = vtkNewPoints->InsertNextPoint(tmp);
container->GetPointIds()->InsertNextId(id);
}
vtkNewCells->InsertNextCell(container);
}
vtkSmartPointer<vtkPolyData> fiberPolyData = vtkSmartPointer<vtkPolyData>::New();
fiberPolyData->SetPoints(vtkNewPoints);
fiberPolyData->SetLines(vtkNewCells);
mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New();
vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New();
matrix->Identity();
if (use_matrix)
for (int i=0; i<4; ++i)
for (int j=0; j<4; ++j)
{
if (j<3)
matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]/m_Header.voxel_size[j]);
else
matrix->SetElement(i, j, m_Header.vox_to_ras[i][j]);
}
if (m_Header.voxel_order[0]=='R')
{
matrix->SetElement(0,0,-matrix->GetElement(0,0));
matrix->SetElement(0,1,-matrix->GetElement(0,1));
matrix->SetElement(0,2,-matrix->GetElement(0,2));
matrix->SetElement(0,3,-matrix->GetElement(0,3));
}
if (m_Header.voxel_order[1]=='A')
{
matrix->SetElement(1,0,-matrix->GetElement(1,0));
matrix->SetElement(1,1,-matrix->GetElement(1,1));
matrix->SetElement(1,2,-matrix->GetElement(1,2));
matrix->SetElement(1,3,-matrix->GetElement(1,3));
}
if (m_Header.voxel_order[2]=='I')
{
matrix->SetElement(2,0,-matrix->GetElement(2,0));
matrix->SetElement(2,1,-matrix->GetElement(2,1));
matrix->SetElement(2,2,-matrix->GetElement(2,2));
matrix->SetElement(2,3,-matrix->GetElement(2,3));
}
+
+
+ if (flip_x)
+ {
+ matrix->SetElement(0,0,-matrix->GetElement(0,0));
+ matrix->SetElement(0,1,-matrix->GetElement(0,1));
+ matrix->SetElement(0,2,-matrix->GetElement(0,2));
+ matrix->SetElement(0,3,-matrix->GetElement(0,3));
+ }
+ if (flip_y)
+ {
+ matrix->SetElement(1,0,-matrix->GetElement(1,0));
+ matrix->SetElement(1,1,-matrix->GetElement(1,1));
+ matrix->SetElement(1,2,-matrix->GetElement(1,2));
+ matrix->SetElement(1,3,-matrix->GetElement(1,3));
+ }
+ if (flip_z)
+ {
+ matrix->SetElement(2,0,-matrix->GetElement(2,0));
+ matrix->SetElement(2,1,-matrix->GetElement(2,1));
+ matrix->SetElement(2,2,-matrix->GetElement(2,2));
+ matrix->SetElement(2,3,-matrix->GetElement(2,3));
+ }
geometry->SetIndexToWorldTransformByVtkMatrix(matrix);
vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
transformFilter->SetInputData(fiberPolyData);
transformFilter->SetTransform(geometry->GetVtkTransform());
transformFilter->Update();
fib->SetFiberPolyData(transformFilter->GetOutput());
fib->SetTrackVisHeader(geometry.GetPointer());
fib->SetTrackVisHeader(m_Header);
return numPoints;
}
// Update the field in the header to the new FIBER TOTAL.
// ------------------------------------------------------
void TrackVisFiberReader::updateTotal( int totFibers )
{
fseek(m_FilePointer, 1000-12, SEEK_SET);
if (fwrite((char*)&totFibers, 1, 4, m_FilePointer) != 4)
MITK_ERROR << "[ERROR] Problems saving the fiber!";
}
void TrackVisFiberReader::writeHdr()
{
fseek(m_FilePointer, 0, SEEK_SET);
if (fwrite((char*)&m_Header, 1, 1000, m_FilePointer) != 1000)
MITK_ERROR << "[ERROR] Problems saving the fiber!";
}
// Close the TrackVis file, but keep the metadata in the header.
// -------------------------------------------------------------
void TrackVisFiberReader::close()
{
fclose(m_FilePointer);
m_FilePointer = nullptr;
}
bool TrackVisFiberReader::IsTransformValid()
{
if (fabs(m_Header.image_orientation_patient[0])<=0.001 || fabs(m_Header.image_orientation_patient[3])<=0.001 || fabs(m_Header.image_orientation_patient[5])<=0.001)
return false;
return true;
}
diff --git a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h
index ceaf409..1a021b6 100644
--- a/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h
+++ b/Modules/DiffusionCore/IODataStructures/mitkTrackvis.h
@@ -1,53 +1,53 @@
/*===================================================================
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 _TRACKVIS
#define _TRACKVIS
#include <mitkCommon.h>
#include <mitkFiberBundle.h>
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkPoints.h>
#include <vtkPolyLine.h>
#include <itkSize.h>
// Class to handle TrackVis files.
// -------------------------------
class MITKDIFFUSIONCORE_EXPORT TrackVisFiberReader
{
private:
std::string m_Filename;
FILE* m_FilePointer;
mitk::FiberBundle::TrackVis_header m_Header;
public:
short create(std::string m_Filename, mitk::FiberBundle *fib, bool print_header);
short open(std::string m_Filename );
- short read(mitk::FiberBundle* fib , bool use_matrix, bool print_header);
+ short read( mitk::FiberBundle* fib, bool use_matrix, bool flip_x, bool flip_y, bool flip_z, bool print_header);
short write(const mitk::FiberBundle* fib );
void writeHdr();
void updateTotal( int totFibers );
void close();
bool IsTransformValid();
void print_header();
TrackVisFiberReader();
~TrackVisFiberReader();
};
#endif
diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp
index ea66f0d..1b22f33 100644
--- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp
+++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp
@@ -1,94 +1,100 @@
/*===================================================================
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 "mitkFiberBundleTrackVisReader.h"
#include <itkMetaDataObject.h>
#include <vtkPolyData.h>
#include <vtkDataReader.h>
#include <vtkPolyDataReader.h>
#include <vtkMatrix4x4.h>
#include <vtkPolyLine.h>
#include <vtkCellArray.h>
#include <vtkDataArray.h>
#include <vtkFloatArray.h>
#include <vtkCellData.h>
#include <vtkPointData.h>
#include <itksys/SystemTools.hxx>
#include <vtkCleanPolyData.h>
#include <mitkTrackvis.h>
#include <mitkCustomMimeType.h>
#include "mitkDiffusionIOMimeTypes.h"
mitk::FiberBundleTrackVisReader::FiberBundleTrackVisReader()
: mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_TRK_MIMETYPE_NAME(), "TrackVis Fiber Bundle Reader" )
{
Options defaultOptions;
defaultOptions["Apply index to world transform stored in the TRK header"] = true;
- defaultOptions["Print header"] = false;
+ defaultOptions["Print header"] = true;
+ defaultOptions["Flip x"] = false;
+ defaultOptions["Flip y"] = false;
+ defaultOptions["Flip z"] = false;
this->SetDefaultOptions(defaultOptions);
m_ServiceReg = this->RegisterService();
}
mitk::FiberBundleTrackVisReader::FiberBundleTrackVisReader(const FiberBundleTrackVisReader &other)
:mitk::AbstractFileReader(other)
{
}
mitk::FiberBundleTrackVisReader * mitk::FiberBundleTrackVisReader::Clone() const
{
return new FiberBundleTrackVisReader(*this);
}
std::vector<itk::SmartPointer<mitk::BaseData> > mitk::FiberBundleTrackVisReader::DoRead()
{
std::vector<itk::SmartPointer<mitk::BaseData> > result;
try
{
const std::string& locale = "C";
const std::string& currLocale = setlocale( LC_ALL, nullptr );
setlocale(LC_ALL, locale.c_str());
std::string filename = this->GetInputLocation();
MITK_INFO << "Loading tractogram (TrackVis format): " << itksys::SystemTools::GetFilenameName(filename);
std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename);
ext = itksys::SystemTools::LowerCase(ext);
if (ext==".trk")
{
Options options = this->GetOptions();
bool apply_matrix = us::any_cast<bool>(options["Apply index to world transform stored in the TRK header"]);
+ bool flip_x = us::any_cast<bool>(options["Flip x"]);
+ bool flip_y = us::any_cast<bool>(options["Flip y"]);
+ bool flip_z = us::any_cast<bool>(options["Flip z"]);
bool print_header = us::any_cast<bool>(options["Print header"]);
FiberBundle::Pointer mitk_fib = FiberBundle::New();
TrackVisFiberReader reader;
reader.open(this->GetInputLocation().c_str());
- reader.read(mitk_fib.GetPointer(), apply_matrix, print_header);
+ reader.read(mitk_fib.GetPointer(), apply_matrix, flip_x, flip_y, flip_z, print_header);
result.push_back(mitk_fib.GetPointer());
return result;
}
setlocale(LC_ALL, currLocale.c_str());
MITK_INFO << "Fiber bundle read";
}
catch(...)
{
throw;
}
return result;
}