diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp index d0ec59a1ac..6a7a366e84 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp @@ -1,291 +1,326 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. 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 "mitkNrrdTensorImageReader.h" #include "itkImageFileReader.h" #include "itkImageRegionIterator.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "mitkITKImageImport.h" #include "mitkImageDataItem.h" namespace mitk { - void NrrdTensorImageReader - ::GenerateData() - { +void NrrdTensorImageReader +::GenerateData() +{ if ( m_FileName == "") { - throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!"); + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!"); } else { - try - { - const std::string& locale = "C"; - const std::string& currLocale = setlocale( LC_ALL, NULL ); - - if ( locale.compare(currLocale)!=0 ) + try { - try - { - setlocale(LC_ALL, locale.c_str()); - } - catch(...) - { - MITK_INFO << "Could not set locale " << locale; - } - } + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, NULL ); - typedef itk::VectorImage ImageType; - itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); - typedef itk::ImageFileReader FileReaderType; - FileReaderType::Pointer reader = FileReaderType::New(); - reader->SetImageIO(io); - reader->SetFileName(this->m_FileName); - reader->Update(); - ImageType::Pointer img = reader->GetOutput(); - - typedef itk::Image,3> VecImgType; - VecImgType::Pointer vecImg = VecImgType::New(); - vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing - vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin - vecImg->SetDirection( img->GetDirection() ); // Set the image direction - vecImg->SetRegions( img->GetLargestPossibleRegion()); - vecImg->Allocate(); - - itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); - ot = ot.Begin(); - - itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); - it = it.Begin(); - - typedef ImageType::PixelType VarPixType; - typedef VecImgType::PixelType FixPixType; - int numComponents = img->GetNumberOfComponentsPerPixel(); - - itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); - std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); - std::vector::const_iterator itKey = imgMetaKeys.begin(); - std::string metaString; - - bool readFrame = false; - double xx, xy, xz, yx, yy, yz, zx, zy, zz; - MeasurementFrameType measFrame; - measFrame.SetIdentity(); - MeasurementFrameType measFrameTransp; - measFrameTransp.SetIdentity(); - - for (; itKey != imgMetaKeys.end(); itKey ++) - { - itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); - if (itKey->find("measurement frame") != std::string::npos) - { - sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); - - if (xx>10e-10 || xy>10e-10 || xz>10e-10 || - yx>10e-10 || yy>10e-10 || yz>10e-10 || - zx>10e-10 || zy>10e-10 || zz>10e-10 ) + if ( locale.compare(currLocale)!=0 ) { - readFrame = true; - - measFrame(0,0) = xx; - measFrame(0,1) = xy; - measFrame(0,2) = xz; - measFrame(1,0) = yx; - measFrame(1,1) = yy; - measFrame(1,2) = yz; - measFrame(2,0) = zx; - measFrame(2,1) = zy; - measFrame(2,2) = zz; - - measFrameTransp = measFrame.GetTranspose(); + try + { + setlocale(LC_ALL, locale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not set locale " << locale; + } } - } - } - - if (numComponents==6) - { - while (!it.IsAtEnd()) - { - // T'=RTR' - VarPixType vec = it.Get(); - FixPixType fixVec(vec.GetDataPointer()); - if(readFrame) + MITK_INFO << "Loading tensor image ..."; + + typedef itk::VectorImage ImageType; + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + typedef itk::ImageFileReader FileReaderType; + FileReaderType::Pointer reader = FileReaderType::New(); + reader->SetImageIO(io); + reader->SetFileName(this->m_FileName); + reader->Update(); + ImageType::Pointer img = reader->GetOutput(); + + typedef itk::Image,3> VecImgType; + VecImgType::Pointer vecImg = VecImgType::New(); + vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing + vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin + vecImg->SetDirection( img->GetDirection() ); // Set the image direction + vecImg->SetRegions( img->GetLargestPossibleRegion()); + vecImg->Allocate(); + + itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); + ot = ot.Begin(); + + itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); + it = it.Begin(); + + typedef ImageType::PixelType VarPixType; + typedef VecImgType::PixelType FixPixType; + int numComponents = img->GetNumberOfComponentsPerPixel(); + + itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); + std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); + std::vector::const_iterator itKey = imgMetaKeys.begin(); + std::string metaString; + + bool readFrame = false; + double xx, xy, xz, yx, yy, yz, zx, zy, zz; + MeasurementFrameType measFrame; + measFrame.SetIdentity(); + MeasurementFrameType measFrameTransp; + measFrameTransp.SetIdentity(); + + for (; itKey != imgMetaKeys.end(); itKey ++) { - itk::DiffusionTensor3D tensor; - tensor.SetElement(0, vec.GetElement(0)); - tensor.SetElement(1, vec.GetElement(1)); - tensor.SetElement(2, vec.GetElement(2)); - tensor.SetElement(3, vec.GetElement(3)); - tensor.SetElement(4, vec.GetElement(4)); - tensor.SetElement(5, vec.GetElement(5)); - - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); - fixVec = tensor; + itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); + if (itKey->find("measurement frame") != std::string::npos) + { + sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); + + if (xx>10e-10 || xy>10e-10 || xz>10e-10 || + yx>10e-10 || yy>10e-10 || yz>10e-10 || + zx>10e-10 || zy>10e-10 || zz>10e-10 ) + { + readFrame = true; + + measFrame(0,0) = xx; + measFrame(0,1) = xy; + measFrame(0,2) = xz; + measFrame(1,0) = yx; + measFrame(1,1) = yy; + measFrame(1,2) = yz; + measFrame(2,0) = zx; + measFrame(2,1) = zy; + measFrame(2,2) = zz; + + measFrameTransp = measFrame.GetTranspose(); + } + } } - ot.Set(fixVec); - ++ot; - ++it; - } - } - else if(numComponents==9) - { - while (!it.IsAtEnd()) - { - VarPixType vec = it.Get(); - itk::DiffusionTensor3D tensor; - tensor.SetElement(0, vec.GetElement(0)); - tensor.SetElement(1, vec.GetElement(1)); - tensor.SetElement(2, vec.GetElement(2)); - tensor.SetElement(3, vec.GetElement(4)); - tensor.SetElement(4, vec.GetElement(5)); - tensor.SetElement(5, vec.GetElement(8)); - - if(readFrame) + if (numComponents==6) + { + while (!it.IsAtEnd()) + { + // T'=RTR' + VarPixType vec = it.Get(); + FixPixType fixVec(vec.GetDataPointer()); + + if(readFrame) + { + itk::DiffusionTensor3D tensor; + tensor.SetElement(0, vec.GetElement(0)); + tensor.SetElement(1, vec.GetElement(1)); + tensor.SetElement(2, vec.GetElement(2)); + tensor.SetElement(3, vec.GetElement(3)); + tensor.SetElement(4, vec.GetElement(4)); + tensor.SetElement(5, vec.GetElement(5)); + + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + fixVec = tensor; + } + + ot.Set(fixVec); + ++ot; + ++it; + } + } + else if(numComponents==9) + { + while (!it.IsAtEnd()) + { + VarPixType vec = it.Get(); + itk::DiffusionTensor3D tensor; + tensor.SetElement(0, vec.GetElement(0)); + tensor.SetElement(1, vec.GetElement(1)); + tensor.SetElement(2, vec.GetElement(2)); + tensor.SetElement(3, vec.GetElement(4)); + tensor.SetElement(4, vec.GetElement(5)); + tensor.SetElement(5, vec.GetElement(8)); + + if(readFrame) + { + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + } + + FixPixType fixVec(tensor); + ot.Set(fixVec); + ++ot; + ++it; + } + } + else if (numComponents==1) { - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); - tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + typedef itk::Image ImageType; + itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); + typedef itk::ImageFileReader FileReaderType; + FileReaderType::Pointer reader = FileReaderType::New(); + reader->SetImageIO(io); + reader->SetFileName(this->m_FileName); + reader->Update(); + ImageType::Pointer img = reader->GetOutput(); + + itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); + + while (!ot.IsAtEnd()) + { + itk::DiffusionTensor3D tensor; + ImageType::IndexType idx; + idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; + for (int te=0; teGetPixel(idx)); + } + if(readFrame) + { + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); + tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); + } + FixPixType fixVec(tensor); + ot.Set(fixVec); + ++ot; + } + } + else + { + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); } - FixPixType fixVec(tensor); - ot.Set(fixVec); - ++ot; - ++it; - } - } - else - { - throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); - } + this->GetOutput()->InitializeByItk(vecImg.GetPointer()); + this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); - this->GetOutput()->InitializeByItk(vecImg.GetPointer()); - this->GetOutput()->SetVolume(vecImg->GetBufferPointer()); + try + { + setlocale(LC_ALL, currLocale.c_str()); + } + catch(...) + { + MITK_INFO << "Could not reset locale " << currLocale; + } - try + } + catch(std::exception& e) { - setlocale(LC_ALL, currLocale.c_str()); + throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { - MITK_INFO << "Could not reset locale " << currLocale; + throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested DTI file!"); } - - } - catch(std::exception& e) - { - throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); - } - catch(...) - { - throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested DTI file!"); - } } - } +} - void NrrdTensorImageReader::GenerateOutputInformation() - { +void NrrdTensorImageReader::GenerateOutputInformation() +{ - } +} - const char* NrrdTensorImageReader - ::GetFileName() const - { +const char* NrrdTensorImageReader +::GetFileName() const +{ return m_FileName.c_str(); - } +} - void NrrdTensorImageReader - ::SetFileName(const char* aFileName) - { +void NrrdTensorImageReader +::SetFileName(const char* aFileName) +{ m_FileName = aFileName; - } +} - const char* NrrdTensorImageReader - ::GetFilePrefix() const - { +const char* NrrdTensorImageReader +::GetFilePrefix() const +{ return m_FilePrefix.c_str(); - } +} - void NrrdTensorImageReader - ::SetFilePrefix(const char* aFilePrefix) - { +void NrrdTensorImageReader +::SetFilePrefix(const char* aFilePrefix) +{ m_FilePrefix = aFilePrefix; - } +} - const char* NrrdTensorImageReader - ::GetFilePattern() const - { +const char* NrrdTensorImageReader +::GetFilePattern() const +{ return m_FilePattern.c_str(); - } +} - void NrrdTensorImageReader - ::SetFilePattern(const char* aFilePattern) - { +void NrrdTensorImageReader +::SetFilePattern(const char* aFilePattern) +{ m_FilePattern = aFilePattern; - } +} - bool NrrdTensorImageReader - ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) - { +bool NrrdTensorImageReader +::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) +{ // First check the extension if( filename == "" ) { - return false; + return false; } std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext == ".hdti" || ext == ".dti") { - return true; + return true; } return false; - } +} - itk::DiffusionTensor3D NrrdTensorImageReader - ::ConvertMatrixTypeToFixedArrayType(const itk::DiffusionTensor3D::Superclass::MatrixType & matrix) - { - /* | 0 1 2 | +itk::DiffusionTensor3D NrrdTensorImageReader +::ConvertMatrixTypeToFixedArrayType(const itk::DiffusionTensor3D::Superclass::MatrixType & matrix) +{ + /* | 0 1 2 | * | X 3 4 | * | X X 5 | */ itk::DiffusionTensor3D arr; arr.SetElement(0,matrix(0,0)); arr.SetElement(1,matrix(0,1)); arr.SetElement(2,matrix(0,2)); arr.SetElement(3,matrix(1,3)); arr.SetElement(4,matrix(1,4)); arr.SetElement(5,matrix(2,5)); return arr; - } +} } //namespace MITK diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.cpp new file mode 100755 index 0000000000..68543c2b6e --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.cpp @@ -0,0 +1,351 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +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 __itkEvaluateTractogramDirectionsFilter_cpp +#define __itkEvaluateTractogramDirectionsFilter_cpp + +#include "itkEvaluateTractogramDirectionsFilter.h" +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +namespace itk { + +template< class PixelType > +EvaluateTractogramDirectionsFilter< PixelType > +::EvaluateTractogramDirectionsFilter(): + m_ReferenceImageSet(NULL), + m_IgnoreMissingDirections(false), + m_Eps(0.0001), + m_UseInterpolation(false) +{ + this->SetNumberOfOutputs(1); +} + +template< class PixelType > +itk::Vector EvaluateTractogramDirectionsFilter< PixelType >::GetItkVector(double point[3]) +{ + itk::Vector itkVector; + itkVector[0] = point[0]; + itkVector[1] = point[1]; + itkVector[2] = point[2]; + return itkVector; +} + +template< class PixelType > +vnl_vector_fixed EvaluateTractogramDirectionsFilter< PixelType >::GetVnlVector(double point[3]) +{ + vnl_vector_fixed vnlVector; + vnlVector[0] = point[0]; + vnlVector[1] = point[1]; + vnlVector[2] = point[2]; + return vnlVector; +} + +template< class PixelType > +vnl_vector_fixed EvaluateTractogramDirectionsFilter< PixelType >::GetVnlVector(Vector& vector) +{ + vnl_vector_fixed vnlVector; + vnlVector[0] = vector[0]; + vnlVector[1] = vector[1]; + vnlVector[2] = vector[2]; + return vnlVector; +} + +template< class PixelType > +itk::Point EvaluateTractogramDirectionsFilter< PixelType >::GetItkPoint(double point[3]) +{ + itk::Point itkPoint; + itkPoint[0] = point[0]; + itkPoint[1] = point[1]; + itkPoint[2] = point[2]; + return itkPoint; +} + +template< class PixelType > +void EvaluateTractogramDirectionsFilter< PixelType >::GenerateData() +{ + if (m_Tractogram.IsNull() || m_ReferenceImageSet.IsNull()) + return; + + if (m_UseInterpolation) + MITK_INFO << "Using trilinear interpolation"; + else + MITK_INFO << "Using nearest neighbor interpolation"; + + if (m_IgnoreMissingDirections) + MITK_INFO << "Ignoring missing directions"; + else + MITK_INFO << "Penalizing missing directions"; + + // angular error image + typename OutputImageType::Pointer outputImage = OutputImageType::New(); + outputImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + outputImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + outputImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + outputImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + outputImage->Allocate(); + outputImage->FillBuffer(0.0); + + DoubleImageType::Pointer counterImage = DoubleImageType::New(); + counterImage->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + counterImage->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + counterImage->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + counterImage->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + counterImage->Allocate(); + counterImage->FillBuffer(0.0); + + std::vector< DoubleImageType::Pointer > directionFound; + for (int i=0; iSize(); i++) + { + DoubleImageType::Pointer check = DoubleImageType::New(); + check->SetOrigin( m_ReferenceImageSet->GetElement(0)->GetOrigin() ); + check->SetRegions( m_ReferenceImageSet->GetElement(0)->GetLargestPossibleRegion() ); + check->SetSpacing( m_ReferenceImageSet->GetElement(0)->GetSpacing() ); + check->SetDirection( m_ReferenceImageSet->GetElement(0)->GetDirection() ); + check->Allocate(); + check->FillBuffer(90); + directionFound.push_back(check); + } + + if (m_MaskImage.IsNull()) + { + m_MaskImage = UCharImageType::New(); + m_MaskImage->SetOrigin( outputImage->GetOrigin() ); + m_MaskImage->SetRegions( outputImage->GetLargestPossibleRegion() ); + m_MaskImage->SetSpacing( outputImage->GetSpacing() ); + m_MaskImage->SetDirection( outputImage->GetDirection() ); + m_MaskImage->Allocate(); + m_MaskImage->FillBuffer(1); + } + + m_MeanAngularError = 0.0; + m_MedianAngularError = 0; + m_MaxAngularError = 0.0; + m_MinAngularError = itk::NumericTraits::max(); + m_VarAngularError = 0.0; + m_AngularErrorVector.clear(); + + float minSpacing = 1; + if(outputImage->GetSpacing()[0]GetSpacing()[1] && outputImage->GetSpacing()[0]GetSpacing()[2]) + minSpacing = outputImage->GetSpacing()[0]; + else if (outputImage->GetSpacing()[1] < outputImage->GetSpacing()[2]) + minSpacing = outputImage->GetSpacing()[1]; + else + minSpacing = outputImage->GetSpacing()[2]; + FiberBundleType::Pointer fiberBundle = m_Tractogram->GetDeepCopy(); + fiberBundle->ResampleFibers(minSpacing/10); + + MITK_INFO << "Evaluating tractogram"; + vtkSmartPointer fiberPolyData = fiberBundle->GetFiberPolyData(); + boost::progress_display disp( fiberBundle->GetNumFibers() ); + for( int i=0; iGetNumFibers(); i++ ) + { + vtkCell* cell = fiberPolyData->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + if (numPoints<2) + continue; + + for( int j=0; jGetPoint(j); + itk::Point vertex = GetItkPoint(temp); + itk::Vector v = GetItkVector(temp); + + itk::Vector dir(3); + if (jGetPoint(j+1))-v; + else + dir = v-GetItkVector(points->GetPoint(j-1)); + + vnl_vector_fixed< PixelType, 3 > fiberDir = GetVnlVector(dir); + fiberDir.normalize(); + + itk::Index<3> idx; + itk::ContinuousIndex contIndex; + m_MaskImage->TransformPhysicalPointToIndex(vertex, idx); + m_MaskImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); + + if (!m_UseInterpolation) // use nearest neighbour interpolation + { + if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)<=0) + continue; + + double angle = 90; + int usedIndex = -1; + for (int k=0; kSize(); k++) // and each test direction + { + vnl_vector_fixed< PixelType, 3 > refDir = m_ReferenceImageSet->GetElement(k)->GetPixel(idx).GetVnlVector(); + if (refDir.magnitude()>m_Eps) // normalize if not null + refDir.normalize(); + else + continue; + + // calculate angle between directions + double tempAngle = acos(fabs(dot_product(refDir, fiberDir)))*180.0/M_PI; + directionFound.at(k)->SetPixel(idx, tempAngle); + + if (tempAngle < angle) + { + angle = tempAngle; + usedIndex = k; + } + } + if (usedIndex>=0) + directionFound.at(usedIndex)->SetPixel(idx, -1); + else if (m_IgnoreMissingDirections) + angle = 0; + + counterImage->SetPixel(idx, counterImage->GetPixel(idx)+1); + outputImage->SetPixel(idx, outputImage->GetPixel(idx)+angle); + continue; + } + + double frac_x = contIndex[0] - idx[0]; + double frac_y = contIndex[1] - idx[1]; + double frac_z = contIndex[2] - idx[2]; + if (frac_x<0) + { + idx[0] -= 1; + frac_x += 1; + } + if (frac_y<0) + { + idx[1] -= 1; + frac_y += 1; + } + if (frac_z<0) + { + idx[2] -= 1; + frac_z += 1; + } + + // use trilinear interpolation + itk::Index<3> newIdx; + for (int x=0; x<2; x++) + { + frac_x = 1-frac_x; + for (int y=0; y<2; y++) + { + frac_y = 1-frac_y; + for (int z=0; z<2; z++) + { + frac_z = 1-frac_z; + + newIdx[0] = idx[0]+x; + newIdx[1] = idx[1]+y; + newIdx[2] = idx[2]+z; + + double frac = frac_x*frac_y*frac_z; + + // is position valid? + if (!m_MaskImage->GetLargestPossibleRegion().IsInside(newIdx) || m_MaskImage->GetPixel(newIdx)<=0) + continue; + + double angle = 90; + int usedIndex = -1; + for (int k=0; kSize(); k++) // and each test direction + { + vnl_vector_fixed< PixelType, 3 > refDir = m_ReferenceImageSet->GetElement(k)->GetPixel(newIdx).GetVnlVector(); + if (refDir.magnitude()>m_Eps) // normalize if not null + refDir.normalize(); + else + continue; + + // calculate angle between directions + double tempAngle = acos(fabs(dot_product(refDir, fiberDir)))*180.0/M_PI; + directionFound.at(k)->SetPixel(newIdx, tempAngle); + + if (tempAngle < angle) + { + angle = tempAngle; + usedIndex = k; + } + } + if (usedIndex>=0) + directionFound.at(usedIndex)->SetPixel(newIdx, -1); + else if (m_IgnoreMissingDirections) + angle = 0; + + counterImage->SetPixel(newIdx, counterImage->GetPixel(newIdx)+1); + outputImage->SetPixel(newIdx, outputImage->GetPixel(newIdx)+frac*angle); + } + } + } + } + ++disp; + } + + ImageRegionIterator< OutputImageType > oit(outputImage, outputImage->GetLargestPossibleRegion()); + ImageRegionIterator< UCharImageType > mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion()); + ImageRegionIterator< DoubleImageType > cit(counterImage, counterImage->GetLargestPossibleRegion()); + + + if (!m_IgnoreMissingDirections) + MITK_INFO << "Creatings statistics and accounting for missing directions"; + else + MITK_INFO << "Creatings statistics"; + boost::progress_display disp2(outputImage->GetLargestPossibleRegion().GetNumberOfPixels()); + while( !oit.IsAtEnd() ) + { + if ( mit.Get()>0) + { + if ( cit.Get()>m_Eps ) + oit.Set(oit.Get()/cit.Get()); + + if (!m_IgnoreMissingDirections) + { + int missingCount = 0; + + if ( cit.Get()>m_Eps ) + missingCount++; + + for (int i=0; iGetPixel(oit.GetIndex())>0 && m_ReferenceImageSet->GetElement(i)->GetPixel(oit.GetIndex()).GetVnlVector().magnitude()>m_Eps ) + { + oit.Set(oit.Get()+directionFound.at(i)->GetPixel(oit.GetIndex())); + missingCount++; + } + } + if (missingCount>1) + oit.Set(oit.Get()/missingCount); + } + + if (oit.Get()>m_MaxAngularError) + m_MaxAngularError = oit.Get(); + if (oit.Get()SetNthOutput(0, outputImage); +} + +} + +#endif // __itkEvaluateTractogramDirectionsFilter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.h new file mode 100755 index 0000000000..cf1cddc6c7 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkEvaluateTractogramDirectionsFilter.h @@ -0,0 +1,107 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +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. + +===================================================================*/ + +/*=================================================================== + +This file is based heavily on a corresponding ITK filter. + +===================================================================*/ +#ifndef __itkEvaluateTractogramDirectionsFilter_h_ +#define __itkEvaluateTractogramDirectionsFilter_h_ + +#include +#include +#include +#include + +namespace itk{ +/** \class EvaluateTractogramDirectionsFilter + */ + +template< class PixelType > +class EvaluateTractogramDirectionsFilter : public ImageSource< Image< PixelType, 3 > > +{ + +public: + + typedef EvaluateTractogramDirectionsFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageSource< Image< PixelType, 3 > > Superclass; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef typename Superclass::OutputImageType OutputImageType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(EvaluateTractogramDirectionsFilter, ImageToImageFilter) + + typedef mitk::FiberBundleX FiberBundleType; + typedef Vector< float, 3 > DirectionType; + typedef Image< DirectionType, 3 > DirectionImageType; + typedef VectorContainer< int, DirectionImageType::Pointer > DirectionImageContainerType; + typedef Image< float, 3 > FloatImageType; + typedef Image< bool, 3 > BoolImageType; + typedef Image< unsigned char, 3 > UCharImageType; + typedef Image< double, 3 > DoubleImageType; + + itkSetMacro( Tractogram, FiberBundleType::Pointer) + itkSetMacro( ReferenceImageSet , DirectionImageContainerType::Pointer) + itkSetMacro( MaskImage , UCharImageType::Pointer) + itkSetMacro( IgnoreMissingDirections , bool) + itkSetMacro( UseInterpolation , bool) + + itkGetMacro( MeanAngularError, float) + itkGetMacro( MinAngularError, float) + itkGetMacro( MaxAngularError, float) + itkGetMacro( VarAngularError, float) + itkGetMacro( MedianAngularError, float) + +protected: + EvaluateTractogramDirectionsFilter(); + ~EvaluateTractogramDirectionsFilter() {} + + void GenerateData(); + + itk::Point GetItkPoint(double point[3]); + itk::Vector GetItkVector(double point[3]); + vnl_vector_fixed GetVnlVector(double point[3]); + vnl_vector_fixed GetVnlVector(Vector< PixelType, 3 >& vector); + + UCharImageType::Pointer m_MaskImage; + DirectionImageContainerType::Pointer m_ReferenceImageSet; + bool m_IgnoreMissingDirections; + double m_MeanAngularError; + double m_MedianAngularError; + double m_MaxAngularError; + double m_MinAngularError; + double m_VarAngularError; + std::vector< double > m_AngularErrorVector; + + double m_Eps; + FiberBundleType::Pointer m_Tractogram; + bool m_UseInterpolation; +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkEvaluateTractogramDirectionsFilter.cpp" +#endif + +#endif //__itkEvaluateTractogramDirectionsFilter_h_ + diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp index 39376f9aea..6907da7f64 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.cpp @@ -1,824 +1,833 @@ #include "itkTractsToVectorImageFilter.h" // VTK #include #include #include // ITK #include #include // misc #define _USE_MATH_DEFINES #include #include namespace itk{ static bool CompareVectorLengths(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { return (v1.magnitude()>v2.magnitude()); } -TractsToVectorImageFilter::TractsToVectorImageFilter(): +template< class PixelType > +TractsToVectorImageFilter< PixelType >::TractsToVectorImageFilter(): m_AngularThreshold(0.7), m_MaskImage(NULL), m_NumDirectionsImage(NULL), m_NormalizeVectors(false), m_Epsilon(0.999), m_UseWorkingCopy(true), m_MaxNumDirections(3), m_UseTrilinearInterpolation(false), m_Thres(0.5) { this->SetNumberOfRequiredOutputs(1); } -TractsToVectorImageFilter::~TractsToVectorImageFilter() +template< class PixelType > +TractsToVectorImageFilter< PixelType >::~TractsToVectorImageFilter() { } -vnl_vector_fixed TractsToVectorImageFilter::GetVnlVector(double point[3]) +template< class PixelType > +vnl_vector_fixed TractsToVectorImageFilter< PixelType >::GetVnlVector(double point[3]) { vnl_vector_fixed vnlVector; vnlVector[0] = point[0]; vnlVector[1] = point[1]; vnlVector[2] = point[2]; return vnlVector; } -itk::Point TractsToVectorImageFilter::GetItkPoint(double point[3]) +template< class PixelType > +itk::Point TractsToVectorImageFilter< PixelType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } -void TractsToVectorImageFilter::GenerateData() +template< class PixelType > +void TractsToVectorImageFilter< PixelType >::GenerateData() { mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); // calculate new image parameters itk::Vector spacing; itk::Point origin; itk::Matrix direction; ImageRegion<3> imageRegion; if (!m_MaskImage.IsNull()) { spacing = m_MaskImage->GetSpacing(); imageRegion = m_MaskImage->GetLargestPossibleRegion(); origin = m_MaskImage->GetOrigin(); direction = m_MaskImage->GetDirection(); } else { spacing = geometry->GetSpacing(); origin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); origin[0] += bounds.GetElement(0); origin[1] += bounds.GetElement(2); origin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) direction[j][i] = geometry->GetMatrixColumn(i)[j]; imageRegion.SetSize(0, geometry->GetExtent(0)); imageRegion.SetSize(1, geometry->GetExtent(1)); imageRegion.SetSize(2, geometry->GetExtent(2)); m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } OutputImageType::RegionType::SizeType outImageSize = imageRegion.GetSize(); m_OutImageSpacing = m_MaskImage->GetSpacing(); m_ClusteredDirectionsContainer = ContainerType::New(); // initialize crossings image m_CrossingsImage = ItkUcharImgType::New(); m_CrossingsImage->SetSpacing( spacing ); m_CrossingsImage->SetOrigin( origin ); m_CrossingsImage->SetDirection( direction ); m_CrossingsImage->SetRegions( imageRegion ); m_CrossingsImage->Allocate(); m_CrossingsImage->FillBuffer(0); // initialize num directions image m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); // resample fiber bundle float minSpacing = 1; if(m_OutImageSpacing[0]GetDeepCopy(); // resample fiber bundle for sufficient voxel coverage m_FiberBundle->ResampleFibers(minSpacing/10); // iterate over all fibers vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numFibers = m_FiberBundle->GetNumFibers(); itk::TimeProbe clock; m_DirectionsContainer = ContainerType::New(); if (m_UseTrilinearInterpolation) MITK_INFO << "Generating directions from tractogram (trilinear interpolation)"; else MITK_INFO << "Generating directions from tractogram"; boost::progress_display disp(numFibers); for( int i=0; iGetNextCell ( numPoints, points ); if (numPoints<2) continue; itk::Index<3> index; index.Fill(0); itk::ContinuousIndex contIndex; vnl_vector_fixed dir, wDir; itk::Point worldPos; vnl_vector v; for( int j=0; jGetPoint(points[j]); worldPos = GetItkPoint(temp); v = GetVnlVector(temp); dir = GetVnlVector(fiberPolyData->GetPoint(points[j+1]))-v; dir.normalize(); m_MaskImage->TransformPhysicalPointToIndex(worldPos, index); m_MaskImage->TransformPhysicalPointToContinuousIndex(worldPos, contIndex); if (m_MaskImage->GetPixel(index)==0) continue; if (!m_UseTrilinearInterpolation) { if (index[0] < 0 || index[0] >= outImageSize[0]) continue; if (index[1] < 0 || index[1] >= outImageSize[1]) continue; if (index[2] < 0 || index[2] >= outImageSize[2]) continue; int idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]); DirectionContainerType::Pointer dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, dir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), dir); } else { dirCont->InsertElement(0, dir); m_DirectionsContainer->InsertElement(idx, dirCont); } continue; } float frac_x = contIndex[0] - index[0]; float frac_y = contIndex[1] - index[1]; float frac_z = contIndex[2] - index[2]; if (frac_x<0) { index[0] -= 1; frac_x += 1; } if (frac_y<0) { index[1] -= 1; frac_y += 1; } if (frac_z<0) { index[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (index[0] < 0 || index[0] >= outImageSize[0]-1) continue; if (index[1] < 0 || index[1] >= outImageSize[1]-1) continue; if (index[2] < 0 || index[2] >= outImageSize[2]-1) continue; DirectionContainerType::Pointer dirCont; int idx; wDir = dir; float weight = ( frac_x)*( frac_y)*( frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2] ); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = ( frac_x)*(1-frac_y)*( frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0] + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2] ); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = ( frac_x)*( frac_y)*(1-frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0] + outImageSize[0]*(index[1] + outImageSize[1]*index[2]+outImageSize[1]); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = ( frac_x)*(1-frac_y)*(1-frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0] + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = (1-frac_x)*( frac_y)*( frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0]+1 + outImageSize[0]*(index[1] + outImageSize[1]*index[2] ); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = (1-frac_x)*( frac_y)*(1-frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0]+1 + outImageSize[0]*(index[1] + outImageSize[1]*index[2]+outImageSize[1]); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = (1-frac_x)*(1-frac_y)*( frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2] ); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } wDir = dir; weight = (1-frac_x)*(1-frac_y)*(1-frac_z); if (weight>m_Thres) { wDir *= weight; idx = index[0]+1 + outImageSize[0]*(index[1]+1+ outImageSize[1]*index[2]+outImageSize[1]); dirCont = DirectionContainerType::New(); if (m_DirectionsContainer->IndexExists(idx)) { dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull()) { dirCont = DirectionContainerType::New(); dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } else dirCont->InsertElement(dirCont->Size(), wDir); } else { dirCont->InsertElement(0, wDir); m_DirectionsContainer->InsertElement(idx, dirCont); } } } clock.Stop(); } vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); itk::ImageRegionIterator dirIt(m_NumDirectionsImage, m_NumDirectionsImage->GetLargestPossibleRegion()); itk::ImageRegionIterator crossIt(m_CrossingsImage, m_CrossingsImage->GetLargestPossibleRegion()); m_DirectionImageContainer = DirectionImageContainerType::New(); int maxNumDirections = 0; MITK_INFO << "Clustering directions"; boost::progress_display disp2(outImageSize[0]*outImageSize[1]*outImageSize[2]); for(crossIt.GoToBegin(); !crossIt.IsAtEnd(); ++crossIt) { ++disp2; OutputImageType::IndexType index = crossIt.GetIndex(); int idx = index[0]+(index[1]+index[2]*outImageSize[1])*outImageSize[0]; if (!m_DirectionsContainer->IndexExists(idx)) { ++dirIt; continue; } DirectionContainerType::Pointer dirCont = m_DirectionsContainer->GetElement(idx); if (dirCont.IsNull() || index[0] < 0 || index[0] >= outImageSize[0] || index[1] < 0 || index[1] >= outImageSize[1] || index[2] < 0 || index[2] >= outImageSize[2]) { ++dirIt; continue; } std::vector< DirectionType > directions; for (int i=0; iSize(); i++) if (dirCont->ElementAt(i).magnitude()>0.0001) directions.push_back(dirCont->ElementAt(i)); if (!directions.empty()) directions = FastClustering(directions); std::sort( directions.begin(), directions.end(), CompareVectorLengths ); if ( directions.size() > maxNumDirections ) { for (int i=maxNumDirections; iSetSpacing( spacing ); directionImage->SetOrigin( origin ); directionImage->SetDirection( direction ); directionImage->SetRegions( imageRegion ); directionImage->Allocate(); Vector< float, 3 > nullVec; nullVec.Fill(0.0); directionImage->FillBuffer(nullVec); m_DirectionImageContainer->InsertElement(i, directionImage); } maxNumDirections = std::min((int)directions.size(), m_MaxNumDirections); } int numDir = directions.size(); if (numDir>m_MaxNumDirections) numDir = m_MaxNumDirections; for (int i=0; i container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); DirectionType dir = m_MaskImage->GetDirection()*directions.at(i); // set direction image pixel ItkDirectionImageType::Pointer directionImage = m_DirectionImageContainer->GetElement(i); Vector< float, 3 > pixel; pixel.SetElement(0, dir[0]); pixel.SetElement(1, dir[1]); pixel.SetElement(2, dir[2]); directionImage->SetPixel(index, pixel); // add direction to vector field (with spacing compensation) itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2*minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2*minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2*minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2*minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2*minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2*minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } dirIt.Set(numDir); ++dirIt; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundleX::New(directionsPolyData); } -std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter::FastClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs) +template< class PixelType > +std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::FastClustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs) { std::vector< vnl_vector_fixed< double, 3 > > outDirs; if (inDirs.empty()) return outDirs; vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean; std::vector< vnl_vector_fixed< double, 3 > > normalizedDirs; std::vector< int > touched; for (int i=0; i0.0001) { counter = 0; oldMean = currentMean; workingMean = oldMean; workingMean.normalize(); currentMean.fill(0.0); for (int i=0; i=m_AngularThreshold) { currentMean += inDirs[i]; touched[i] = 1; counter++; } else if (-angle>=m_AngularThreshold) { currentMean -= inDirs[i]; touched[i] = 1; counter++; } } } // found stable mean if (counter>0) { currentMean /= counter; float mag = currentMean.magnitude(); if (mag>0) { if (mag>max) max = mag; outDirs.push_back(currentMean); } } // find next unused seed free = false; for (int i=0; i0) for (int i=0; i > TractsToVectorImageFilter::Clustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs) +template< class PixelType > +std::vector< vnl_vector_fixed< double, 3 > > TractsToVectorImageFilter< PixelType >::Clustering(std::vector< vnl_vector_fixed< double, 3 > >& inDirs) { std::vector< vnl_vector_fixed< double, 3 > > outDirs; if (inDirs.empty()) return outDirs; vnl_vector_fixed< double, 3 > oldMean, currentMean, workingMean; std::vector< vnl_vector_fixed< double, 3 > > normalizedDirs; std::vector< int > touched; for (int i=0; i0.0001) { counter = 0; oldMean = currentMean; workingMean = oldMean; workingMean.normalize(); currentMean.fill(0.0); for (int i=0; i=m_AngularThreshold) { currentMean += inDirs[i]; counter++; } else if (-angle>=m_AngularThreshold) { currentMean -= inDirs[i]; counter++; } } } // found stable mean if (counter>0) { bool add = true; vnl_vector_fixed< double, 3 > normMean = currentMean; normMean.normalize(); for (int i=0; i dir = outDirs[i]; dir.normalize(); if ((normMean-dir).magnitude()<=0.0001) { add = false; break; } } currentMean /= counter; if (add) { float mag = currentMean.magnitude(); if (mag>0) { if (mag>max) max = mag; outDirs.push_back(currentMean); } } } } if (m_NormalizeVectors) for (int i=0; i0) for (int i=0; i +TractsToVectorImageFilter< PixelType >::DirectionContainerType::Pointer TractsToVectorImageFilter< PixelType >::MeanShiftClustering(DirectionContainerType::Pointer dirCont) { DirectionContainerType::Pointer container = DirectionContainerType::New(); float max = 0; for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) { vnl_vector_fixed mean = ClusterStep(dirCont, it.Value()); if (mean.is_zero()) continue; bool addMean = true; for (DirectionContainerType::ConstIterator it2 = container->Begin(); it2!=container->End(); ++it2) { vnl_vector_fixed dir = it2.Value(); float angle = fabs(dot_product(mean, dir)/(mean.magnitude()*dir.magnitude())); if (angle>=m_Epsilon) { addMean = false; break; } } if (addMean) { if (m_NormalizeVectors) mean.normalize(); else if (mean.magnitude()>max) max = mean.magnitude(); container->InsertElement(container->Size(), mean); } } // max normalize voxel directions if (max>0 && !m_NormalizeVectors) for (int i=0; iSize(); i++) container->ElementAt(i) /= max; if (container->Size()Size()) return MeanShiftClustering(container); else return container; } -vnl_vector_fixed TractsToVectorImageFilter::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed currentMean) +template< class PixelType > +vnl_vector_fixed TractsToVectorImageFilter< PixelType >::ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed currentMean) { vnl_vector_fixed newMean; newMean.fill(0); for (DirectionContainerType::ConstIterator it = dirCont->Begin(); it!=dirCont->End(); ++it) { vnl_vector_fixed dir = it.Value(); float angle = dot_product(currentMean, dir)/(currentMean.magnitude()*dir.magnitude()); if (angle>=m_AngularThreshold) newMean += dir; else if (-angle>=m_AngularThreshold) newMean -= dir; } if (fabs(dot_product(currentMean, newMean)/(currentMean.magnitude()*newMean.magnitude()))>=m_Epsilon || newMean.is_zero()) return newMean; else return ClusterStep(dirCont, newMean); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h index 0b8746cfe6..f7fd1e434e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToVectorImageFilter.h @@ -1,107 +1,108 @@ #ifndef __itkTractsToVectorImageFilter_h__ #define __itkTractsToVectorImageFilter_h__ // MITK #include // ITK #include #include // VTK #include #include #include #include #include using namespace mitk; namespace itk{ /** * \brief Extracts the voxel-wise main directions of the input fiber bundle. */ +template< class PixelType > class TractsToVectorImageFilter : public ImageSource< VectorImage< float, 3 > > { public: typedef TractsToVectorImageFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef itk::Vector OutputVectorType; - typedef itk::Image OutputImageType; + typedef itk::Image OutputImageType; typedef std::vector< OutputImageType::Pointer > OutputImageContainerType; typedef vnl_vector_fixed< double, 3 > DirectionType; typedef VectorContainer< int, DirectionType > DirectionContainerType; typedef VectorContainer< int, DirectionContainerType::Pointer > ContainerType; typedef Image< Vector< float, 3 >, 3> ItkDirectionImageType; typedef VectorContainer< int, ItkDirectionImageType::Pointer > DirectionImageContainerType; typedef itk::Image ItkUcharImgType; itkNewMacro(Self) itkTypeMacro( TractsToVectorImageFilter, ImageSource ) itkSetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold itkGetMacro( AngularThreshold, float) ///< cluster directions that are closer together than the specified threshold itkSetMacro( NormalizeVectors, bool) ///< Normalize vectors to length 1 itkGetMacro( NormalizeVectors, bool) ///< Normalize vectors to length 1 itkSetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. itkGetMacro( UseWorkingCopy, bool) ///< Do not modify input fiber bundle. Use a copy. itkSetMacro( MaxNumDirections, int) ///< If more directions are extracted, only the largest are kept. itkGetMacro( MaxNumDirections, int) ///< If more directions are extracted, only the largest are kept. itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< only process voxels inside mask itkSetMacro( FiberBundle, FiberBundleX::Pointer) ///< input fiber bundle itkGetMacro( ClusteredDirectionsContainer, ContainerType::Pointer) ///< output directions itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< nimber of directions per voxel itkGetMacro( CrossingsImage, ItkUcharImgType::Pointer) ///< mask voxels containing crossings itkGetMacro( OutputFiberBundle, FiberBundleX::Pointer) ///< vector field for visualization purposes itkGetMacro( DirectionImageContainer, DirectionImageContainerType::Pointer) ///< output directions void GenerateData(); protected: std::vector< DirectionType > Clustering(std::vector< DirectionType >& inDirs); std::vector< DirectionType > FastClustering(std::vector< DirectionType >& inDirs); ///< cluster fiber directions DirectionContainerType::Pointer MeanShiftClustering(DirectionContainerType::Pointer dirCont); vnl_vector_fixed ClusterStep(DirectionContainerType::Pointer dirCont, vnl_vector_fixed currentMean); vnl_vector_fixed GetVnlVector(double point[3]); itk::Point GetItkPoint(double point[3]); TractsToVectorImageFilter(); virtual ~TractsToVectorImageFilter(); FiberBundleX::Pointer m_FiberBundle; ///< input fiber bundle float m_AngularThreshold; ///< cluster directions that are closer together than the specified threshold float m_Epsilon; ///< epsilon for vector equality check ItkUcharImgType::Pointer m_MaskImage; ///< only voxels inside the binary mask are processed bool m_NormalizeVectors; ///< normalize vectors to length 1 itk::Vector m_OutImageSpacing; ///< spacing of output image ContainerType::Pointer m_DirectionsContainer; ///< container for fiber directions bool m_UseWorkingCopy; ///< do not modify input fiber bundle but work on copy bool m_UseTrilinearInterpolation; ///< trilinearly interpolate between neighbouring voxels int m_MaxNumDirections; ///< if more directions per voxel are extracted, only the largest are kept float m_Thres; ///< distance threshold for trilinear interpolation // output datastructures ContainerType::Pointer m_ClusteredDirectionsContainer; ///< contains direction vectors for each voxel ItkUcharImgType::Pointer m_NumDirectionsImage; ///< shows number of fibers per voxel ItkUcharImgType::Pointer m_CrossingsImage; ///< shows voxels containing more than one fiber DirectionImageContainerType::Pointer m_DirectionImageContainer; ///< contains images that contain the output directions FiberBundleX::Pointer m_OutputFiberBundle; ///< vector field for visualization purposes }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractsToVectorImageFilter.cpp" #endif #endif // __itkTractsToVectorImageFilter_h__ diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt index d8ab51b688..0828ba2dd0 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/CMakeLists.txt @@ -1,51 +1,54 @@ OPTION(BUILD_FiberTrackingMiniApps "Build commandline tools for fiber tracking" OFF) IF(BUILD_FiberTrackingMiniApps) # include necessary modules here MitkExt QmitkExt MITK_CHECK_MODULE(_RESULT DiffusionCore FiberTracking ) IF(_RESULT) MESSAGE("Warning: FiberTrackingMiniApps is missing ${_RESULT}") ELSE(_RESULT) MITK_USE_MODULE( DiffusionCore FiberTracking ) # needed include directories INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${ALL_INCLUDE_DIRECTORIES}) PROJECT( mitkFiberTrackingMiniApps ) # fill in the standalone executables here SET(FIBERTRACKINGMINIAPPS mitkFiberTrackingMiniApps ) # set additional files here SET(FIBERTRACKING_ADDITIONAL_FILES MiniAppManager.cpp GibbsTracking.cpp StreamlineTracking.cpp FiberProcessing.cpp PeakExtraction.cpp TractometerAngularErrorTool.cpp + TractogramAngularError.cpp + FiberDirectionExtraction.cpp + ImageFormatConverter.cpp ) # create an executable foreach tool (only one at the moment) FOREACH(tool ${FIBERTRACKINGMINIAPPS}) ADD_EXECUTABLE( ${tool} ${tool}.cpp ${FIBERTRACKING_ADDITIONAL_FILES} ) TARGET_LINK_LIBRARIES( ${tool} ${ALL_LIBRARIES} ) ENDFOREACH(tool) ENDIF() MITK_INSTALL_TARGETS(EXECUTABLES mitkFiberTrackingMiniApps ) ENDIF(BUILD_FiberTrackingMiniApps) diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/FiberDirectionExtraction.cpp similarity index 51% copy from Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp copy to Modules/DiffusionImaging/FiberTracking/MiniApps/FiberDirectionExtraction.cpp index b7f099d30e..3186c7c108 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/FiberDirectionExtraction.cpp @@ -1,264 +1,167 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. 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 "MiniAppManager.h" #include #include #include #include #include #include -#include #include #include "ctkCommandLineParser.h" -#include "ctkCommandLineParser.cpp" #include #include -#include #include #include -#include -#include +#include +#include #define _USE_MATH_DEFINES #include -int TractometerAngularErrorTool(int argc, char* argv[]) +int FiberDirectionExtraction(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", mitk::Any(), false); - parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", mitk::Any(), false); parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", mitk::Any(), false); parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); - parser.addArgument("ignore", "n", ctkCommandLineParser::Bool, "don't increase error for missing or too many directions"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; - ctkCommandLineParser::StringContainerType referenceImages = mitk::any_cast(parsedArgs["reference"]); - string fibFile = mitk::any_cast(parsedArgs["input"]); string maskImage(""); if (parsedArgs.count("mask")) maskImage = mitk::any_cast(parsedArgs["mask"]); float angularThreshold = 25; if (parsedArgs.count("athresh")) angularThreshold = mitk::any_cast(parsedArgs["athresh"]); string outRoot = mitk::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = mitk::any_cast(parsedArgs["verbose"]); - bool ignore = false; - if (parsedArgs.count("ignore")) - ignore = mitk::any_cast(parsedArgs["ignore"]); try { RegisterDiffusionCoreObjectFactory(); RegisterFiberTrackingObjectFactory(); typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; - typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); - // load reference directions - ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); - for (int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); - typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; - CasterType::Pointer caster = CasterType::New(); - caster->SetInput(img); - caster->Update(); - ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); - referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); - } - catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } - } - // load/create mask image - ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); - if (maskImage.compare("")==0) - { - ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); - itkMaskImage->SetSpacing( dirImg->GetSpacing() ); - itkMaskImage->SetOrigin( dirImg->GetOrigin() ); - itkMaskImage->SetDirection( dirImg->GetDirection() ); - itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); - itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); - itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); - itkMaskImage->Allocate(); - itkMaskImage->FillBuffer(1); - } - else + ItkUcharImgType::Pointer itkMaskImage = NULL; + if (maskImage.compare("")!=0) { + MITK_INFO << "Using mask image"; + itkMaskImage = ItkUcharImgType::New(); mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); - CastToItkImage(mitkMaskImage, itkMaskImage); + mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } - // extract directions from fiber bundle - itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); fOdfFilter->SetNormalizeVectors(true); fOdfFilter->SetUseWorkingCopy(false); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); + // write direction images + for (int i=0; iSize(); i++) + { + itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); + typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; + WriterType::Pointer writer = WriterType::New(); + + string outfilename = outRoot; + outfilename.append("_DIRECTION_"); + outfilename.append(boost::lexical_cast(i)); + outfilename.append(".nrrd"); + + MITK_INFO << "writing " << outfilename; + writer->SetFileName(outfilename.c_str()); + writer->SetInput(itkImg); + writer->Update(); + } + if (verbose) { // write vector field mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle(); mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) { string outfilename = outRoot; outfilename.append("_VECTOR_FIELD.fib"); (*it)->SetFileName( outfilename.c_str() ); (*it)->DoWrite( directions.GetPointer() ); } } - // write direction images - for (int i=0; iSize(); i++) - { - itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); - typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; - WriterType::Pointer writer = WriterType::New(); - - string outfilename = outRoot; - outfilename.append("_DIRECTION_"); - outfilename.append(boost::lexical_cast(i)); - outfilename.append(".nrrd"); - - MITK_INFO << "writing " << outfilename; - writer->SetFileName(outfilename.c_str()); - writer->SetInput(itkImg); - writer->Update(); - } - // write num direction image { ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS.nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } - // evaluate directions - EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); - evaluationFilter->SetImageSet(directionImageContainer); - evaluationFilter->SetReferenceImageSet(referenceImageContainer); - evaluationFilter->SetMaskImage(itkMaskImage); - evaluationFilter->SetIgnoreMissingDirections(ignore); - evaluationFilter->Update(); - - if (verbose) - { - EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); - typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; - WriterType::Pointer writer = WriterType::New(); - - string outfilename = outRoot; - outfilename.append("_ERROR_IMAGE.nrrd"); - - MITK_INFO << "writing " << outfilename; - writer->SetFileName(outfilename.c_str()); - writer->SetInput(angularErrorImage); - writer->Update(); - } - - string logFile = outRoot; - logFile.append("_ANGULAR_ERROR.csv"); - - ofstream file; - file.open (logFile.c_str()); - - string sens = "Mean:"; - sens.append(","); - sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); - sens.append(";\n"); - - sens.append("Median:"); - sens.append(","); - sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); - sens.append(";\n"); - - sens.append("Maximum:"); - sens.append(","); - sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); - sens.append(";\n"); - - sens.append("Minimum:"); - sens.append(","); - sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); - sens.append(";\n"); - - sens.append("STDEV:"); - sens.append(","); - sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); - sens.append(";\n"); - - file << sens; - - file.close(); - MITK_INFO << "DONE"; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } -RegisterFiberTrackingMiniApp(TractometerAngularErrorTool); +RegisterFiberTrackingMiniApp(FiberDirectionExtraction); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/ImageFormatConverter.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/ImageFormatConverter.cpp new file mode 100755 index 0000000000..0b0d1fb97e --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/ImageFormatConverter.cpp @@ -0,0 +1,84 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +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 "MiniAppManager.h" +#include +#include +#include +#include +#include +#include + +using namespace mitk; +#include "ctkCommandLineParser.h" + +int ImageFormatConverter(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("in", "i", ctkCommandLineParser::String, "input image", mitk::Any(), false); + parser.addArgument("out", "o", ctkCommandLineParser::String, "output image", mitk::Any(), false); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + // mandatory arguments + string imageName = mitk::any_cast(parsedArgs["in"]); + string outImage = mitk::any_cast(parsedArgs["out"]); + + try + { + RegisterDiffusionCoreObjectFactory(); + + MITK_INFO << "Loading image ..."; + const std::string s1="", s2=""; + std::vector infile = BaseDataIO::LoadBaseDataFromFile( imageName, s1, s2, false ); + + MITK_INFO << "Writing " << outImage; + if ( dynamic_cast*>(infile.at(0).GetPointer()) ) + { + DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); + NrrdDiffusionImageWriter::Pointer writer = NrrdDiffusionImageWriter::New(); + writer->SetFileName(outImage); + writer->SetInput(dwi); + writer->Update(); + } + else if ( dynamic_cast(infile.at(0).GetPointer()) ) + { + Image::Pointer image = dynamic_cast(infile.at(0).GetPointer()); + mitk::IOUtil::SaveImage(image, outImage); + } + } + catch (itk::ExceptionObject e) + { + MITK_INFO << e; + return EXIT_FAILURE; + } + catch (std::exception e) + { + MITK_INFO << e.what(); + return EXIT_FAILURE; + } + catch (...) + { + MITK_INFO << "ERROR!?!"; + return EXIT_FAILURE; + } + MITK_INFO << "DONE"; + return EXIT_SUCCESS; +} +RegisterFiberTrackingMiniApp(ImageFormatConverter); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractogramAngularError.cpp similarity index 67% copy from Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp copy to Modules/DiffusionImaging/FiberTracking/MiniApps/TractogramAngularError.cpp index b7f099d30e..5e70d879c7 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractogramAngularError.cpp @@ -1,264 +1,205 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. 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 "MiniAppManager.h" #include #include #include #include #include #include -#include #include #include "ctkCommandLineParser.h" -#include "ctkCommandLineParser.cpp" #include #include -#include #include #include #include #include +#include #define _USE_MATH_DEFINES #include -int TractometerAngularErrorTool(int argc, char* argv[]) +int TractogramAngularError(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", mitk::Any(), false); parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", mitk::Any(), false); parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", mitk::Any(), false); parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); - parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); parser.addArgument("ignore", "n", ctkCommandLineParser::Bool, "don't increase error for missing or too many directions"); + parser.addArgument("trilinear", "t", ctkCommandLineParser::Bool, "use trilinear instead of nearest neighbor interpolation"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; ctkCommandLineParser::StringContainerType referenceImages = mitk::any_cast(parsedArgs["reference"]); string fibFile = mitk::any_cast(parsedArgs["input"]); string maskImage(""); if (parsedArgs.count("mask")) maskImage = mitk::any_cast(parsedArgs["mask"]); - float angularThreshold = 25; - if (parsedArgs.count("athresh")) - angularThreshold = mitk::any_cast(parsedArgs["athresh"]); - string outRoot = mitk::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = mitk::any_cast(parsedArgs["verbose"]); bool ignore = false; if (parsedArgs.count("ignore")) ignore = mitk::any_cast(parsedArgs["ignore"]); + bool interpolate = false; + if (parsedArgs.count("interpolate")) + interpolate = mitk::any_cast(parsedArgs["interpolate"]); + try { RegisterDiffusionCoreObjectFactory(); RegisterFiberTrackingObjectFactory(); typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; - typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; + typedef itk::EvaluateTractogramDirectionsFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); + if (!inputTractogram) + return EXIT_FAILURE; // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } } // load/create mask image ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); if (maskImage.compare("")==0) { ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); } else { mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); - CastToItkImage(mitkMaskImage, itkMaskImage); - } - - - // extract directions from fiber bundle - itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); - fOdfFilter->SetFiberBundle(inputTractogram); - fOdfFilter->SetMaskImage(itkMaskImage); - fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); - fOdfFilter->SetNormalizeVectors(true); - fOdfFilter->SetUseWorkingCopy(false); - fOdfFilter->Update(); - ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); - - if (verbose) - { - // write vector field - mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle(); - mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); - for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) - { - if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) { - string outfilename = outRoot; - outfilename.append("_VECTOR_FIELD.fib"); - (*it)->SetFileName( outfilename.c_str() ); - (*it)->DoWrite( directions.GetPointer() ); - } - } - - // write direction images - for (int i=0; iSize(); i++) - { - itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); - typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; - WriterType::Pointer writer = WriterType::New(); - - string outfilename = outRoot; - outfilename.append("_DIRECTION_"); - outfilename.append(boost::lexical_cast(i)); - outfilename.append(".nrrd"); - - MITK_INFO << "writing " << outfilename; - writer->SetFileName(outfilename.c_str()); - writer->SetInput(itkImg); - writer->Update(); - } - - // write num direction image - { - ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); - typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; - WriterType::Pointer writer = WriterType::New(); - - string outfilename = outRoot; - outfilename.append("_NUM_DIRECTIONS.nrrd"); - - MITK_INFO << "writing " << outfilename; - writer->SetFileName(outfilename.c_str()); - writer->SetInput(numDirImage); - writer->Update(); - } + mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); - evaluationFilter->SetImageSet(directionImageContainer); + evaluationFilter->SetTractogram(inputTractogram); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); evaluationFilter->SetIgnoreMissingDirections(ignore); + evaluationFilter->SetUseInterpolation(interpolate); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string logFile = outRoot; logFile.append("_ANGULAR_ERROR.csv"); ofstream file; file.open (logFile.c_str()); string sens = "Mean:"; sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(";\n"); sens.append("Median:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(";\n"); sens.append("Maximum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(";\n"); sens.append("Minimum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(";\n"); sens.append("STDEV:"); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; file.close(); MITK_INFO << "DONE"; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } -RegisterFiberTrackingMiniApp(TractometerAngularErrorTool); +RegisterFiberTrackingMiniApp(TractogramAngularError); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp index b7f099d30e..78297e6f0c 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/TractometerAngularErrorTool.cpp @@ -1,264 +1,264 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. 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 "MiniAppManager.h" #include #include #include #include #include #include #include #include #include "ctkCommandLineParser.h" #include "ctkCommandLineParser.cpp" +#include #include #include -#include #include #include #include #include #define _USE_MATH_DEFINES #include int TractometerAngularErrorTool(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", ctkCommandLineParser::String, "input tractogram (.fib, vtk ascii file format)", mitk::Any(), false); parser.addArgument("reference", "r", ctkCommandLineParser::StringList, "reference direction images", mitk::Any(), false); parser.addArgument("out", "o", ctkCommandLineParser::String, "output root", mitk::Any(), false); parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); parser.addArgument("athresh", "a", ctkCommandLineParser::Float, "angular threshold in degrees. closer fiber directions are regarded as one direction and clustered together.", 25, true); parser.addArgument("verbose", "v", ctkCommandLineParser::Bool, "output optional and intermediate calculation results"); parser.addArgument("ignore", "n", ctkCommandLineParser::Bool, "don't increase error for missing or too many directions"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; ctkCommandLineParser::StringContainerType referenceImages = mitk::any_cast(parsedArgs["reference"]); string fibFile = mitk::any_cast(parsedArgs["input"]); string maskImage(""); if (parsedArgs.count("mask")) maskImage = mitk::any_cast(parsedArgs["mask"]); float angularThreshold = 25; if (parsedArgs.count("athresh")) angularThreshold = mitk::any_cast(parsedArgs["athresh"]); string outRoot = mitk::any_cast(parsedArgs["out"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = mitk::any_cast(parsedArgs["verbose"]); bool ignore = false; if (parsedArgs.count("ignore")) ignore = mitk::any_cast(parsedArgs["ignore"]); try { RegisterDiffusionCoreObjectFactory(); RegisterFiberTrackingObjectFactory(); typedef itk::Image ItkUcharImgType; typedef itk::Image< itk::Vector< float, 3>, 3 > ItkDirectionImage3DType; typedef itk::VectorContainer< int, ItkDirectionImage3DType::Pointer > ItkDirectionImageContainerType; typedef itk::EvaluateDirectionImagesFilter< float > EvaluationFilterType; // load fiber bundle mitk::FiberBundleX::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::LoadDataNode(fibFile)->GetData()); // load reference directions ItkDirectionImageContainerType::Pointer referenceImageContainer = ItkDirectionImageContainerType::New(); for (int i=0; i(mitk::IOUtil::LoadDataNode(referenceImages.at(i))->GetData()); typedef mitk::ImageToItk< ItkDirectionImage3DType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(img); caster->Update(); ItkDirectionImage3DType::Pointer itkImg = caster->GetOutput(); referenceImageContainer->InsertElement(referenceImageContainer->Size(),itkImg); } catch(...){ MITK_INFO << "could not load: " << referenceImages.at(i); } } // load/create mask image ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); if (maskImage.compare("")==0) { ItkDirectionImage3DType::Pointer dirImg = referenceImageContainer->GetElement(0); itkMaskImage->SetSpacing( dirImg->GetSpacing() ); itkMaskImage->SetOrigin( dirImg->GetOrigin() ); itkMaskImage->SetDirection( dirImg->GetDirection() ); itkMaskImage->SetLargestPossibleRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetBufferedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->SetRequestedRegion( dirImg->GetLargestPossibleRegion() ); itkMaskImage->Allocate(); itkMaskImage->FillBuffer(1); } else { mitk::Image::Pointer mitkMaskImage = dynamic_cast(mitk::IOUtil::LoadDataNode(maskImage)->GetData()); - CastToItkImage(mitkMaskImage, itkMaskImage); + mitk::CastToItkImage(mitkMaskImage, itkMaskImage); } // extract directions from fiber bundle - itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); + itk::TractsToVectorImageFilter::Pointer fOdfFilter = itk::TractsToVectorImageFilter::New(); fOdfFilter->SetFiberBundle(inputTractogram); fOdfFilter->SetMaskImage(itkMaskImage); fOdfFilter->SetAngularThreshold(cos(angularThreshold*M_PI/180)); fOdfFilter->SetNormalizeVectors(true); fOdfFilter->SetUseWorkingCopy(false); fOdfFilter->Update(); ItkDirectionImageContainerType::Pointer directionImageContainer = fOdfFilter->GetDirectionImageContainer(); if (verbose) { // write vector field mitk::FiberBundleX::Pointer directions = fOdfFilter->GetOutputFiberBundle(); mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(directions.GetPointer()) ) { string outfilename = outRoot; outfilename.append("_VECTOR_FIELD.fib"); (*it)->SetFileName( outfilename.c_str() ); (*it)->DoWrite( directions.GetPointer() ); } } // write direction images for (int i=0; iSize(); i++) { - itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); - typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; + itk::TractsToVectorImageFilter::ItkDirectionImageType::Pointer itkImg = directionImageContainer->GetElement(i); + typedef itk::ImageFileWriter< itk::TractsToVectorImageFilter::ItkDirectionImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(itkImg); writer->Update(); } // write num direction image { ItkUcharImgType::Pointer numDirImage = fOdfFilter->GetNumDirectionsImage(); typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_NUM_DIRECTIONS.nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(numDirImage); writer->Update(); } } // evaluate directions EvaluationFilterType::Pointer evaluationFilter = EvaluationFilterType::New(); evaluationFilter->SetImageSet(directionImageContainer); evaluationFilter->SetReferenceImageSet(referenceImageContainer); evaluationFilter->SetMaskImage(itkMaskImage); evaluationFilter->SetIgnoreMissingDirections(ignore); evaluationFilter->Update(); if (verbose) { EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; WriterType::Pointer writer = WriterType::New(); string outfilename = outRoot; outfilename.append("_ERROR_IMAGE.nrrd"); MITK_INFO << "writing " << outfilename; writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } string logFile = outRoot; logFile.append("_ANGULAR_ERROR.csv"); ofstream file; file.open (logFile.c_str()); string sens = "Mean:"; sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(";\n"); sens.append("Median:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(";\n"); sens.append("Maximum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(";\n"); sens.append("Minimum:"); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(";\n"); sens.append("STDEV:"); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append(";\n"); file << sens; file.close(); MITK_INFO << "DONE"; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } RegisterFiberTrackingMiniApp(TractometerAngularErrorTool); diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 641f614cee..ab9a5b3c39 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,105 +1,106 @@ set(CPP_FILES MiniApps/ctkCommandLineParser.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp # IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures IODataStructures/mitkFiberTrackingObjectFactory.cpp # Rendering Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp # Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp #Rendering/mitkPlanarFigureMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp ) set(H_FILES MiniApps/ctkCommandLineParser.h # Rendering Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h # Rendering/mitkFiberBundleXThreadMonitorMapper3D.h #Rendering/mitkPlanarFigureMapper3D.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h # IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h IODataStructures/mitkFiberTrackingObjectFactory.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h Algorithms/itkFibersFromPlanarFiguresFilter.h Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkKspaceImageFilter.h Algorithms/itkDftImageFilter.h Algorithms/itkAddArtifactsToDwiImageFilter.h Algorithms/itkFieldmapGeneratorFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h + Algorithms/itkEvaluateTractogramDirectionsFilter.h # (old) Tractography Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/itkStreamlineTrackingFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h # Signal Models SignalModels/mitkDiffusionSignalModel.h SignalModels/mitkTensorModel.h SignalModels/mitkBallModel.h SignalModels/mitkDotModel.h SignalModels/mitkAstroStickModel.h SignalModels/mitkStickModel.h SignalModels/mitkDiffusionNoiseModel.h SignalModels/mitkRicianNoiseModel.h SignalModels/mitkKspaceArtifact.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml )