diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/MiniApps/CMakeLists.txt index 40fc96fe31..569d92b2c5 100755 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/CMakeLists.txt @@ -1,50 +1,51 @@ OPTION(BUILD_DiffusionCoreMiniApps "Build commandline tools for diffusion core" OFF) IF(BUILD_DiffusionCoreMiniApps) # include necessary modules here MitkExt QmitkExt MITK_CHECK_MODULE(_RESULT DiffusionCore ) IF(_RESULT) MESSAGE("Warning: DiffusionCoreMiniApps is missing ${_RESULT}") ELSE(_RESULT) MITK_USE_MODULE( DiffusionCore ) # needed include directories INCLUDE_DIRECTORIES( ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${ALL_INCLUDE_DIRECTORIES}) PROJECT( mitkDiffusionCoreMiniApps ) # fill in the standalone executables here SET(DIFFUSIONCOREMINIAPPS mitkDiffusionCoreMiniApps ) # set additional files here SET(DIFFUSIONCORE_ADDITIONAL_FILES MiniAppManager.cpp ImageFormatConverter.cpp TensorReconstruction.cpp QballReconstruction.cpp DiffusionIndices.cpp + CopyGeometry.cpp ) # create an executable foreach tool (only one at the moment) FOREACH(tool ${DIFFUSIONCOREMINIAPPS}) ADD_EXECUTABLE( ${tool} ${tool}.cpp ${DIFFUSIONCORE_ADDITIONAL_FILES} ) TARGET_LINK_LIBRARIES( ${tool} ${ALL_LIBRARIES} ) ENDFOREACH(tool) ENDIF() MITK_INSTALL_TARGETS(EXECUTABLES mitkDiffusionCoreMiniApps ) ENDIF(BUILD_DiffusionCoreMiniApps) diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/CopyGeometry.cpp similarity index 67% copy from Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp copy to Modules/DiffusionImaging/DiffusionCore/MiniApps/CopyGeometry.cpp index 9e5ee24fc4..35b9ef5e89 100755 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/CopyGeometry.cpp @@ -1,86 +1,93 @@ /*=================================================================== 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 "ctkCommandLineParser.h" -#include "ctkCommandLineParser.cpp" using namespace mitk; #include "ctkCommandLineParser.h" -int ImageFormatConverter(int argc, char* argv[]) +int CopyGeometry(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("in", "i", ctkCommandLineParser::String, "input image", us::Any(), false); + parser.addArgument("ref", "r", ctkCommandLineParser::String, "reference image", us::Any(), false); parser.addArgument("out", "o", ctkCommandLineParser::String, "output image", us::Any(), false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["in"]); + string refImage = us::any_cast(parsedArgs["ref"]); string outImage = us::any_cast(parsedArgs["out"]); try { RegisterDiffusionCoreObjectFactory(); - MITK_INFO << "Loading image ..."; + MITK_INFO << "Loading image " << imageName; const std::string s1="", s2=""; - std::vector infile = BaseDataIO::LoadBaseDataFromFile( imageName, s1, s2, false ); + std::vector infile = BaseDataIO::LoadBaseDataFromFile( refImage, s1, s2, false ); + Image::Pointer source = dynamic_cast(infile.at(0).GetPointer()); + infile = BaseDataIO::LoadBaseDataFromFile( imageName, s1, s2, false ); + Image::Pointer target = dynamic_cast(infile.at(0).GetPointer()); - MITK_INFO << "Writing " << outImage; - if ( dynamic_cast*>(infile.at(0).GetPointer()) ) + mitk::Geometry3D* s_geom = source->GetGeometry(); + mitk::Geometry3D* t_geom = target->GetGeometry(); + + t_geom->SetIndexToWorldTransform(s_geom->GetIndexToWorldTransform()); + target->SetGeometry(t_geom); + + if ( dynamic_cast*>(target.GetPointer()) ) { - DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); + MITK_INFO << "Writing " << outImage; + DiffusionImage::Pointer dwi = dynamic_cast*>(target.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); - } + else + mitk::IOUtil::SaveImage(target, 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; } -RegisterDiffusionCoreMiniApp(ImageFormatConverter); +RegisterDiffusionCoreMiniApp(CopyGeometry); diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/DiffusionIndices.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/DiffusionIndices.cpp index 99b9b5ee58..c01941e9fc 100644 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/DiffusionIndices.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/DiffusionIndices.cpp @@ -1,137 +1,141 @@ /*=================================================================== 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 #include "ctkCommandLineParser.h" #include /** * Calculate indices derived from Qball or tensor images */ int DiffusionIndices(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", ctkCommandLineParser::String, "input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image)", us::Any(), false); parser.addArgument("index", "idx", ctkCommandLineParser::String, "index (fa, gfa, ra, ad, rd, ca, l2, l3, md)", us::Any(), false); parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output file", us::Any(), false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string inFileName = us::any_cast(parsedArgs["input"]); string index = us::any_cast(parsedArgs["index"]); string outFileName = us::any_cast(parsedArgs["outFile"]); try { RegisterDiffusionCoreObjectFactory(); // load input image const std::string s1="", s2=""; std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( inFileName, s1, s2, false ); if( boost::algorithm::ends_with(inFileName, ".qbi") && index=="gfa" ) { typedef itk::Vector OdfVectorType; typedef itk::Image ItkQballImageType; mitk::QBallImage::Pointer mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); ItkQballImageType::Pointer itk_qbi = ItkQballImageType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itk_qbi); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); + MITK_INFO << "Writing " << outFileName; + itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(gfaFilter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if( boost::algorithm::ends_with(inFileName, ".dti") ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itk_dti.GetPointer() ); if(index=="fa") measurementsCalculator->SetMeasure(MeasurementsType::FA); else if(index=="ra") measurementsCalculator->SetMeasure(MeasurementsType::RA); else if(index=="ad") measurementsCalculator->SetMeasure(MeasurementsType::AD); else if(index=="rd") measurementsCalculator->SetMeasure(MeasurementsType::RD); else if(index=="ca") measurementsCalculator->SetMeasure(MeasurementsType::CA); else if(index=="l2") measurementsCalculator->SetMeasure(MeasurementsType::L2); else if(index=="l3") measurementsCalculator->SetMeasure(MeasurementsType::L3); else if(index=="md") measurementsCalculator->SetMeasure(MeasurementsType::MD); else { MITK_WARN << "No valid diffusion index for input image (tensor image) defined"; return EXIT_FAILURE; } measurementsCalculator->Update(); + MITK_INFO << "Writing " << outFileName; + itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(measurementsCalculator->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } } 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; } RegisterDiffusionCoreMiniApp(DiffusionIndices); diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp index 9e5ee24fc4..7add2929c5 100755 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/ImageFormatConverter.cpp @@ -1,86 +1,86 @@ /*=================================================================== 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 "ctkCommandLineParser.h" #include "ctkCommandLineParser.cpp" using namespace mitk; #include "ctkCommandLineParser.h" int ImageFormatConverter(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("in", "i", ctkCommandLineParser::String, "input image", us::Any(), false); parser.addArgument("out", "o", ctkCommandLineParser::String, "output image", us::Any(), false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["in"]); string outImage = us::any_cast(parsedArgs["out"]); try { RegisterDiffusionCoreObjectFactory(); - MITK_INFO << "Loading image ..."; + MITK_INFO << "Loading image " << imageName; 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()) ) { + MITK_INFO << "Writing " << outImage; 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; } RegisterDiffusionCoreMiniApp(ImageFormatConverter); diff --git a/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp b/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp index 45cadfb53b..f011965583 100644 --- a/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/MiniApps/QballReconstruction.cpp @@ -1,111 +1,176 @@ /*=================================================================== 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 #include #include "itkVectorContainer.h" #include "vnl/vnl_vector_fixed.h" #include "itkVectorImage.h" #include "mitkBaseDataIOFactory.h" #include "mitkDiffusionImage.h" #include "mitkBaseData.h" #include "mitkDiffusionCoreObjectFactory.h" #include "mitkCoreObjectFactory.h" #include "mitkCoreExtObjectFactory.h" #include "mitkImageWriter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include #include using namespace mitk; /** * Convert files from one ending to the other */ int QballReconstruction(int argc, char* argv[]) { - if ( argc!=6 - ) + if ( argc<6 ) { std::cout << std::endl; std::cout << "Perform QBall reconstruction on an dwi file" << std::endl; std::cout << std::endl; - std::cout << "usage: " << argv[0] << " " << argv[1] << " " << std::endl; + std::cout << "usage: " << argv[0] << " " << argv[1] << " " << std::endl; std::cout << std::endl; return EXIT_FAILURE; } try { RegisterDiffusionCoreObjectFactory(); MITK_INFO << "Loading image ..."; const std::string s1="", s2=""; std::vector infile = BaseDataIO::LoadBaseDataFromFile( argv[2], s1, s2, false ); DiffusionImage::Pointer dwi = dynamic_cast*>(infile.at(0).GetPointer()); - - typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; - dwi->AverageRedundantGradients(0.001); - FilterType::Pointer filter = FilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetB_Value()); - filter->SetThreshold( boost::lexical_cast(argv[4]) ); - filter->SetLambda(boost::lexical_cast(argv[3])); - filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); - filter->Update(); - mitk::QBallImage::Pointer image = mitk::QBallImage::New(); - image->InitializeByItk( filter->GetOutput() ); - image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + + + int shOrder = 4; + if ( argc>6 ) + shOrder = boost::lexical_cast(argv[6]); + MITK_INFO << "Using SH order of " << shOrder; + switch ( shOrder ) + { + case 4: + { + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); + filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); + filter->Update(); + image->InitializeByItk( filter->GetOutput() ); + image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + break; + } + case 6: + { + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); + filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); + filter->Update(); + image->InitializeByItk( filter->GetOutput() ); + image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + break; + } + case 8: + { + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); + filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); + filter->Update(); + image->InitializeByItk( filter->GetOutput() ); + image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + break; + } + case 10: + { + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); + filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); + filter->Update(); + image->InitializeByItk( filter->GetOutput() ); + image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + break; + } + default: + { + MITK_INFO << "Supplied SH order not supported. Using default order of 4."; + typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); + filter->SetBValue(dwi->GetB_Value()); + filter->SetThreshold( boost::lexical_cast(argv[4]) ); + filter->SetLambda(boost::lexical_cast(argv[3])); + filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); + filter->Update(); + image->InitializeByItk( filter->GetOutput() ); + image->SetVolume( filter->GetOutput()->GetBufferPointer() ); + } + } std::string outfilename = argv[5]; MITK_INFO << "writing image " << outfilename; mitk::NrrdQBallImageWriter::Pointer writer = mitk::NrrdQBallImageWriter::New(); writer->SetInput(image.GetPointer()); writer->SetFileName(outfilename.c_str()); writer->Update(); } catch ( itk::ExceptionObject &err) { MITK_INFO << "Exception: " << err; } catch ( std::exception err) { MITK_INFO << "Exception: " << err.what(); } catch ( ... ) { MITK_INFO << "Exception!"; } return EXIT_SUCCESS; } RegisterDiffusionCoreMiniApp(QballReconstruction); diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index b6f405949d..f6cc685822 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,720 +1,722 @@ /*=================================================================== 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 __itkStreamlineTrackingFilter_txx #define __itkStreamlineTrackingFilter_txx #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #define _USE_MATH_DEFINES #include namespace itk { //#define QBALL_RECON_PI M_PI template< class TTensorPixelType, class TPDPixelType> StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::StreamlineTrackingFilter(): m_FaThreshold(0.2), m_StepSize(1), m_MaxLength(10000), m_SeedsPerVoxel(1), m_F(1.0), m_G(0.0), m_Interpolate(true), m_MinTractLength(0.0) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TTensorPixelType, class TPDPixelType> double StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::BeforeThreadedGenerateData() { m_FiberPolyData = FiberPolyDataType::New(); m_Points = vtkPoints::New(); m_Cells = vtkCellArray::New(); m_InputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); m_ImageSize.resize(3); m_ImageSize[0] = m_InputImage->GetLargestPossibleRegion().GetSize()[0]; m_ImageSize[1] = m_InputImage->GetLargestPossibleRegion().GetSize()[1]; m_ImageSize[2] = m_InputImage->GetLargestPossibleRegion().GetSize()[2]; if (m_ImageSize[0]<3 || m_ImageSize[1]<3 || m_ImageSize[2]<3) m_Interpolate = false; m_ImageSpacing.resize(3); m_ImageSpacing[0] = m_InputImage->GetSpacing()[0]; m_ImageSpacing[1] = m_InputImage->GetSpacing()[1]; m_ImageSpacing[2] = m_InputImage->GetSpacing()[2]; float minSpacing; if(m_ImageSpacing[0]::New(); for (int i=0; iGetNumberOfThreads(); i++) { FiberPolyDataType poly = FiberPolyDataType::New(); m_PolyDataContainer->InsertElement(i, poly); } if (m_SeedImage.IsNull()) { // initialize mask image m_SeedImage = ItkUcharImgType::New(); m_SeedImage->SetSpacing( m_InputImage->GetSpacing() ); m_SeedImage->SetOrigin( m_InputImage->GetOrigin() ); m_SeedImage->SetDirection( m_InputImage->GetDirection() ); m_SeedImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( m_InputImage->GetSpacing() ); m_MaskImage->SetOrigin( m_InputImage->GetOrigin() ); m_MaskImage->SetDirection( m_InputImage->GetDirection() ); m_MaskImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_InputImage->GetSpacing() ); m_FaImage->SetOrigin( m_InputImage->GetOrigin() ); m_FaImage->SetDirection( m_InputImage->GetDirection() ); m_FaImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); m_PdImage = ItkPDImgType::New(); m_PdImage->SetSpacing( m_InputImage->GetSpacing() ); m_PdImage->SetOrigin( m_InputImage->GetOrigin() ); m_PdImage->SetDirection( m_InputImage->GetDirection() ); m_PdImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_PdImage->Allocate(); m_EmaxImage = ItkFloatImgType::New(); m_EmaxImage->SetSpacing( m_InputImage->GetSpacing() ); m_EmaxImage->SetOrigin( m_InputImage->GetOrigin() ); m_EmaxImage->SetDirection( m_InputImage->GetDirection() ); m_EmaxImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_EmaxImage->Allocate(); m_EmaxImage->FillBuffer(1.0); typedef itk::DiffusionTensor3D TensorType; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; for (int x=0; xGetPixel(index); vnl_vector_fixed dir; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); dir.normalize(); m_PdImage->SetPixel(index, dir); m_FaImage->SetPixel(index, tensor.GetFractionalAnisotropy()); m_EmaxImage->SetPixel(index, 2/eigenvalues[2]); } if (m_Interpolate) std::cout << "StreamlineTrackingFilter: using trilinear interpolation" << std::endl; else { if (m_MinCurvatureRadius<0.0) m_MinCurvatureRadius = 0.1*minSpacing; std::cout << "StreamlineTrackingFilter: using nearest neighbor interpolation" << std::endl; } if (m_MinCurvatureRadius<0.0) m_MinCurvatureRadius = 0.5*minSpacing; + m_PointPistance = 0; + std::cout << "StreamlineTrackingFilter: Min. curvature radius: " << m_MinCurvatureRadius << std::endl; std::cout << "StreamlineTrackingFilter: FA threshold: " << m_FaThreshold << std::endl; std::cout << "StreamlineTrackingFilter: stepsize: " << m_StepSize << " mm" << std::endl; std::cout << "StreamlineTrackingFilter: f: " << m_F << std::endl; std::cout << "StreamlineTrackingFilter: g: " << m_G << std::endl; std::cout << "StreamlineTrackingFilter: starting streamline tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::CalculateNewPosition(itk::ContinuousIndex& pos, vnl_vector_fixed& dir, typename InputImageType::IndexType& index) { vnl_matrix_fixed< double, 3, 3 > rot = m_InputImage->GetDirection().GetTranspose(); dir = rot*dir; if (true) { dir *= m_StepSize; pos[0] += dir[0]/m_ImageSpacing[0]; pos[1] += dir[1]/m_ImageSpacing[1]; pos[2] += dir[2]/m_ImageSpacing[2]; index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); } else { dir[0] /= m_ImageSpacing[0]; dir[1] /= m_ImageSpacing[1]; dir[2] /= m_ImageSpacing[2]; int smallest = 0; float x = 100000; if (dir[0]>0) { if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps) x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0]; else x = fabs(pos[0]-std::ceil(pos[0])-0.5)/dir[0]; } else if (dir[0]<0) { if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps) x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0]; else x = -fabs(pos[0]-std::floor(pos[0])+0.5)/dir[0]; } float s = x; float y = 100000; if (dir[1]>0) { if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps) y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1]; else y = fabs(pos[1]-std::ceil(pos[1])-0.5)/dir[1]; } else if (dir[1]<0) { if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps) y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1]; else y = -fabs(pos[1]-std::floor(pos[1])+0.5)/dir[1]; } if (s>y) { s=y; smallest = 1; } float z = 100000; if (dir[2]>0) { if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps) z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2]; else z = fabs(pos[2]-std::ceil(pos[2])-0.5)/dir[2]; } else if (dir[2]<0) { if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps) z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2]; else z = -fabs(pos[2]-std::floor(pos[2])+0.5)/dir[2]; } if (s>z) { s=z; smallest = 2; } // MITK_INFO << "---------------------------------------------"; // MITK_INFO << "s: " << s; // MITK_INFO << "dir: " << dir; // MITK_INFO << "old: " << pos[0] << ", " << pos[1] << ", " << pos[2]; pos[0] += dir[0]*s; pos[1] += dir[1]*s; pos[2] += dir[2]*s; switch (smallest) { case 0: if (dir[0]<0) index[0] = std::floor(pos[0]); else index[0] = std::ceil(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); break; case 1: if (dir[1]<0) index[1] = std::floor(pos[1]); else index[1] = std::ceil(pos[1]); index[0] = RoundToNearest(pos[0]); index[2] = RoundToNearest(pos[2]); break; case 2: if (dir[2]<0) index[2] = std::floor(pos[2]); else index[2] = std::ceil(pos[2]); index[1] = RoundToNearest(pos[1]); index[0] = RoundToNearest(pos[0]); } // float x = 100000; // if (dir[0]>0) // x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0]; // else if (dir[0]<0) // x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0]; // float s = x; // float y = 100000; // if (dir[1]>0) // y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1]; // else if (dir[1]<0) // y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1]; // if (s>y) // s=y; // float z = 100000; // if (dir[2]>0) // z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2]; // else if (dir[2]<0) // z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2]; // if (s>z) // s=z; // s *= 1.001; // pos[0] += dir[0]*s; // pos[1] += dir[1]*s; // pos[2] += dir[2]*s; // index[0] = RoundToNearest(pos[0]); // index[1] = RoundToNearest(pos[1]); // index[2] = RoundToNearest(pos[2]); // MITK_INFO << "new: " << pos[0] << ", " << pos[1] << ", " << pos[2]; } } template< class TTensorPixelType, class TPDPixelType> bool StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::IsValidPosition(itk::ContinuousIndex& pos, typename InputImageType::IndexType &index, vnl_vector_fixed< float, 8 >& interpWeights) { if (!m_InputImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) return false; if (m_Interpolate) { float frac_x = pos[0] - index[0]; float frac_y = pos[1] - index[1]; float frac_z = pos[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] >= m_ImageSize[0]-1) return false; if (index[1] < 0 || index[1] >= m_ImageSize[1]-1) return false; if (index[2] < 0 || index[2] >= m_ImageSize[2]-1) return false; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); typename InputImageType::IndexType tmpIdx; float FA = m_FaImage->GetPixel(index) * interpWeights[0]; tmpIdx = index; tmpIdx[0]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = index; tmpIdx[1]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = index; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[7]; if (FAGetPixel(index) float StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::FollowStreamline(itk::ContinuousIndex pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids) { float tractLength = 0; typedef itk::DiffusionTensor3D TensorType; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; vnl_vector_fixed< float, 8 > interpWeights; typename InputImageType::IndexType index, indexOld; indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1; itk::Point worldPos; float distance = 0; float distanceInVoxel = 0; // starting index and direction index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); vnl_vector_fixed dir = m_PdImage->GetPixel(index); dir *= dirSign; // reverse direction vnl_vector_fixed dirOld = dir; if (dir.magnitude()TransformContinuousIndexToPhysicalPoint( pos, worldPos ); ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); return tractLength; } else if (distance>=m_PointPistance) { m_InputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); distance = 0; } if (!m_Interpolate) // use nearest neighbour interpolation { if (indexOld!=index) // did we enter a new voxel? if yes, calculate new direction { dir = m_PdImage->GetPixel(index); // get principal direction typename InputImageType::PixelType tensor = m_InputImage->GetPixel(index); float scale = m_EmaxImage->GetPixel(index); dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); dir.normalize(); float angle = dot_product(dirOld, dir); if (angle<0) { dir *= -1; angle *= -1; } float r = m_StepSize/(2*std::asin(std::acos(angle)/2)); if (rGetPixel(index) * interpWeights[0]; typename InputImageType::IndexType tmpIdx = index; tmpIdx[0]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = index; tmpIdx[1]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = index; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[7]; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); dir.normalize(); float scale = 2/eigenvalues[2]; dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); dir.normalize(); float angle = dot_product(dirOld, dir); if (angle<0) { dir *= -1; angle *= -1; } float r = m_StepSize/(2*std::asin(std::acos(angle)/2)); if (r void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId) { FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId); vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer Cells = vtkSmartPointer::New(); typedef itk::DiffusionTensor3D TensorType; typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; typedef ImageRegionConstIterator< ItkFloatImgType > FloatIteratorType; typedef typename InputImageType::PixelType InputTensorType; InputIteratorType it(m_InputImage, outputRegionForThread ); MaskIteratorType mit(m_SeedImage, outputRegionForThread ); FloatIteratorType fit(m_FaImage, outputRegionForThread ); MaskIteratorType mit2(m_MaskImage, outputRegionForThread ); it.GoToBegin(); mit.GoToBegin(); mit2.GoToBegin(); fit.GoToBegin(); itk::Point worldPos; while( !it.IsAtEnd() ) { if (mit.Value()==0 || fit.Value() line = vtkSmartPointer::New(); std::vector< vtkIdType > pointIDs; typename InputImageType::IndexType index = it.GetIndex(); itk::ContinuousIndex start; unsigned int counter = 0; if (m_SeedsPerVoxel>1) { start[0] = index[0]+(double)(rand()%99-49)/100; start[1] = index[1]+(double)(rand()%99-49)/100; start[2] = index[2]+(double)(rand()%99-49)/100; } else { start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; } // forward tracking float tractLength = FollowStreamline(start, 1, points, pointIDs); // add ids to line counter += pointIDs.size(); while (!pointIDs.empty()) { line->GetPointIds()->InsertNextId(pointIDs.back()); pointIDs.pop_back(); } // insert start point m_InputImage->TransformContinuousIndexToPhysicalPoint( start, worldPos ); line->GetPointIds()->InsertNextId(points->InsertNextPoint(worldPos.GetDataPointer())); // backward tracking tractLength += FollowStreamline(start, -1, points, pointIDs); counter += pointIDs.size(); if (tractLengthGetPointIds()->InsertNextId(pointIDs.at(i)); Cells->InsertNextCell(line); } ++mit; ++mit2; ++it; ++fit; } poly->SetPoints(points); poly->SetLines(Cells); std::cout << "Thread " << threadId << " finished tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = poly1->GetLines(); vtkSmartPointer vNewPoints = poly1->GetPoints(); vtkSmartPointer vLines = poly2->GetLines(); vLines->InitTraversal(); for( int i=0; iGetNumberOfCells(); i++ ) { vtkIdType numPoints(0); vtkIdType* points(NULL); vLines->GetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(poly2->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); return vNewPolyData; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AfterThreadedGenerateData() { MITK_INFO << "Generating polydata "; m_FiberPolyData = m_PolyDataContainer->GetElement(0); for (int i=1; iGetNumberOfThreads(); i++) { m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i)); } MITK_INFO << "done"; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::PrintSelf(std::ostream& os, Indent indent) const { } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp index 9128879c63..1736e3d8a8 100755 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractsToDWIImageFilter.cpp @@ -1,738 +1,743 @@ /*=================================================================== 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 "itkTractsToDWIImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk { template< class PixelType > TractsToDWIImageFilter< PixelType >::TractsToDWIImageFilter() : m_CircleDummy(false) , m_VolumeAccuracy(10) , m_Upsampling(1) , m_NumberOfRepetitions(1) , m_EnforcePureFiberVoxels(false) , m_InterpolationShrink(1000) , m_FiberRadius(0) , m_SignalScale(25) , m_kOffset(0) , m_tLine(1) , m_UseInterpolation(false) , m_SimulateRelaxation(true) , m_tInhom(50) , m_TE(100) , m_FrequencyMap(NULL) , m_EddyGradientStrength(0.001) , m_SimulateEddyCurrents(false) , m_Spikes(0) , m_SpikeAmplitude(1) { m_Spacing.Fill(2.5); m_Origin.Fill(0.0); m_DirectionMatrix.SetIdentity(); m_ImageRegion.SetSize(0, 10); m_ImageRegion.SetSize(1, 10); m_ImageRegion.SetSize(2, 10); } template< class PixelType > TractsToDWIImageFilter< PixelType >::~TractsToDWIImageFilter() { } template< class PixelType > TractsToDWIImageFilter< PixelType >::DoubleDwiType::Pointer TractsToDWIImageFilter< PixelType >::DoKspaceStuff( std::vector< DoubleDwiType::Pointer >& images ) { // create slice object ImageRegion<2> sliceRegion; sliceRegion.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); sliceRegion.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); Vector< double, 2 > sliceSpacing; sliceSpacing[0] = m_UpsampledSpacing[0]; sliceSpacing[1] = m_UpsampledSpacing[1]; // frequency map slice SliceType::Pointer fMap = NULL; if (m_FrequencyMap.IsNotNull()) { fMap = SliceType::New(); ImageRegion<2> region; region.SetSize(0, m_UpsampledImageRegion.GetSize()[0]); region.SetSize(1, m_UpsampledImageRegion.GetSize()[1]); fMap->SetLargestPossibleRegion( region ); fMap->SetBufferedRegion( region ); fMap->SetRequestedRegion( region ); fMap->Allocate(); } DoubleDwiType::Pointer newImage = DoubleDwiType::New(); newImage->SetSpacing( m_Spacing ); newImage->SetOrigin( m_Origin ); newImage->SetDirection( m_DirectionMatrix ); newImage->SetLargestPossibleRegion( m_ImageRegion ); newImage->SetBufferedRegion( m_ImageRegion ); newImage->SetRequestedRegion( m_ImageRegion ); newImage->SetVectorLength( images.at(0)->GetVectorLength() ); newImage->Allocate(); MatrixType transform = m_DirectionMatrix; for (int i=0; i<3; i++) for (int j=0; j<3; j++) { if (j<2) transform[i][j] *= m_UpsampledSpacing[j]; else transform[i][j] *= m_Spacing[j]; } std::vector< int > spikeVolume; for (int i=0; iGetVectorLength()); std::sort (spikeVolume.begin(), spikeVolume.end()); std::reverse (spikeVolume.begin(), spikeVolume.end()); boost::progress_display disp(images.at(0)->GetVectorLength()*images.at(0)->GetLargestPossibleRegion().GetSize(2)); for (int g=0; gGetVectorLength(); g++) { std::vector< int > spikeSlice; while (!spikeVolume.empty() && spikeVolume.back()==g) { spikeSlice.push_back(rand()%images.at(0)->GetLargestPossibleRegion().GetSize(2)); spikeVolume.pop_back(); } std::sort (spikeSlice.begin(), spikeSlice.end()); std::reverse (spikeSlice.begin(), spikeSlice.end()); for (int z=0; zGetLargestPossibleRegion().GetSize(2); z++) { std::vector< SliceType::Pointer > compartmentSlices; std::vector< double > t2Vector; for (int i=0; i* signalModel; if (iSetLargestPossibleRegion( sliceRegion ); slice->SetBufferedRegion( sliceRegion ); slice->SetRequestedRegion( sliceRegion ); slice->SetSpacing(sliceSpacing); slice->Allocate(); slice->FillBuffer(0.0); // extract slice from channel g for (int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; slice->SetPixel(index2D, images.at(i)->GetPixel(index3D)[g]); if (fMap.IsNotNull() && i==0) fMap->SetPixel(index2D, m_FrequencyMap->GetPixel(index3D)); } compartmentSlices.push_back(slice); t2Vector.push_back(signalModel->GetT2()); } // create k-sapce (inverse fourier transform slices) itk::Size<2> outSize; outSize.SetElement(0, m_ImageRegion.GetSize(0)); outSize.SetElement(1, m_ImageRegion.GetSize(1)); itk::KspaceImageFilter< SliceType::PixelType >::Pointer idft = itk::KspaceImageFilter< SliceType::PixelType >::New(); idft->SetCompartmentImages(compartmentSlices); idft->SetT2(t2Vector); idft->SetkOffset(m_kOffset); idft->SettLine(m_tLine); idft->SetTE(m_TE); idft->SetTinhom(m_tInhom); idft->SetSimulateRelaxation(m_SimulateRelaxation); idft->SetSimulateEddyCurrents(m_SimulateEddyCurrents); idft->SetEddyGradientMagnitude(m_EddyGradientStrength); idft->SetZ((double)z-(double)images.at(0)->GetLargestPossibleRegion().GetSize(2)/2.0); idft->SetDirectionMatrix(transform); idft->SetDiffusionGradientDirection(m_FiberModels.at(0)->GetGradientDirection(g)); idft->SetFrequencyMap(fMap); idft->SetSignalScale(m_SignalScale); idft->SetOutSize(outSize); int numSpikes = 0; while (!spikeSlice.empty() && spikeSlice.back()==z) { numSpikes++; spikeSlice.pop_back(); } idft->SetSpikes(numSpikes); idft->SetSpikeAmplitude(m_SpikeAmplitude); idft->Update(); ComplexSliceType::Pointer fSlice; fSlice = idft->GetOutput(); for (int i=0; iAddArtifact(fSlice); // fourier transform slice SliceType::Pointer newSlice; itk::DftImageFilter< SliceType::PixelType >::Pointer dft = itk::DftImageFilter< SliceType::PixelType >::New(); dft->SetInput(fSlice); dft->Update(); newSlice = dft->GetOutput(); // put slice back into channel g for (int y=0; yGetLargestPossibleRegion().GetSize(1); y++) for (int x=0; xGetLargestPossibleRegion().GetSize(0); x++) { DoubleDwiType::IndexType index3D; index3D[0]=x; index3D[1]=y; index3D[2]=z; SliceType::IndexType index2D; index2D[0]=x; index2D[1]=y; DoubleDwiType::PixelType pix3D = newImage->GetPixel(index3D); pix3D[g] = newSlice->GetPixel(index2D); newImage->SetPixel(index3D, pix3D); } ++disp; } } return newImage; } //template< class PixelType > //TractsToDWIImageFilter< PixelType >::ComplexSliceType::Pointer TractsToDWIImageFilter< PixelType >::RearrangeSlice(ComplexSliceType::Pointer slice) //{ // ImageRegion<2> region = slice->GetLargestPossibleRegion(); // ComplexSliceType::Pointer rearrangedSlice = ComplexSliceType::New(); // rearrangedSlice->SetLargestPossibleRegion( region ); // rearrangedSlice->SetBufferedRegion( region ); // rearrangedSlice->SetRequestedRegion( region ); // rearrangedSlice->Allocate(); // int xHalf = region.GetSize(0)/2; // int yHalf = region.GetSize(1)/2; // for (int y=0; y pix = slice->GetPixel(idx); // if( idx[0] < xHalf ) // idx[0] = idx[0] + xHalf; // else // idx[0] = idx[0] - xHalf; // if( idx[1] < yHalf ) // idx[1] = idx[1] + yHalf; // else // idx[1] = idx[1] - yHalf; // rearrangedSlice->SetPixel(idx, pix); // } // return rearrangedSlice; //} template< class PixelType > void TractsToDWIImageFilter< PixelType >::GenerateData() { // check input data if (m_FiberBundle.IsNull()) itkExceptionMacro("Input fiber bundle is NULL!"); int numFibers = m_FiberBundle->GetNumFibers(); if (numFibers<=0) itkExceptionMacro("Input fiber bundle contains no fibers!"); if (m_FiberModels.empty()) itkExceptionMacro("No diffusion model for fiber compartments defined!"); + if (m_EnforcePureFiberVoxels) + while (m_FiberModels.size()>1) + m_FiberModels.pop_back(); + if (m_NonFiberModels.empty()) itkExceptionMacro("No diffusion model for non-fiber compartments defined!"); int baselineIndex = m_FiberModels[0]->GetFirstBaselineIndex(); if (baselineIndex<0) itkExceptionMacro("No baseline index found!"); // check k-space undersampling if (m_Upsampling<1) m_Upsampling = 1; if (m_TissueMask.IsNotNull()) { // use input tissue mask m_Spacing = m_TissueMask->GetSpacing(); m_Origin = m_TissueMask->GetOrigin(); m_DirectionMatrix = m_TissueMask->GetDirection(); m_ImageRegion = m_TissueMask->GetLargestPossibleRegion(); if (m_Upsampling>1.00001) { MITK_INFO << "Adding ringing artifacts. Image upsampling factor: " << m_Upsampling; ImageRegion<3> region = m_ImageRegion; region.SetSize(0, m_ImageRegion.GetSize(0)*m_Upsampling); region.SetSize(1, m_ImageRegion.GetSize(1)*m_Upsampling); itk::Vector spacing = m_Spacing; spacing[0] /= m_Upsampling; spacing[1] /= m_Upsampling; itk::RescaleIntensityImageFilter::Pointer rescaler = itk::RescaleIntensityImageFilter::New(); rescaler->SetInput(0,m_TissueMask); rescaler->SetOutputMaximum(100); rescaler->SetOutputMinimum(0); rescaler->Update(); itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(rescaler->GetOutput()); resampler->SetOutputParametersFromImage(m_TissueMask); resampler->SetSize(region.GetSize()); resampler->SetOutputSpacing(spacing); resampler->Update(); m_TissueMask = resampler->GetOutput(); } MITK_INFO << "Using tissue mask"; } // initialize output dwi image typename OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( m_Spacing ); outImage->SetOrigin( m_Origin ); outImage->SetDirection( m_DirectionMatrix ); outImage->SetLargestPossibleRegion( m_ImageRegion ); outImage->SetBufferedRegion( m_ImageRegion ); outImage->SetRequestedRegion( m_ImageRegion ); outImage->SetVectorLength( m_FiberModels[0]->GetNumGradients() ); outImage->Allocate(); typename OutputImageType::PixelType temp; temp.SetSize(m_FiberModels[0]->GetNumGradients()); temp.Fill(0.0); outImage->FillBuffer(temp); // is input slize size a power of two? int x=m_ImageRegion.GetSize(0); int y=m_ImageRegion.GetSize(1); if ( x%2 == 1 ) x += 1; if ( y%2 == 1 ) y += 1; // if not, adjust size and dimension (needed for FFT); zero-padding if (x!=m_ImageRegion.GetSize(0)) m_ImageRegion.SetSize(0, x); if (y!=m_ImageRegion.GetSize(1)) m_ImageRegion.SetSize(1, y); // apply undersampling to image parameters m_UpsampledSpacing = m_Spacing; m_UpsampledImageRegion = m_ImageRegion; m_UpsampledSpacing[0] /= m_Upsampling; m_UpsampledSpacing[1] /= m_Upsampling; m_UpsampledImageRegion.SetSize(0, m_ImageRegion.GetSize()[0]*m_Upsampling); m_UpsampledImageRegion.SetSize(1, m_ImageRegion.GetSize()[1]*m_Upsampling); // everything from here on is using the upsampled image parameters!!! if (m_TissueMask.IsNull()) { m_TissueMask = ItkUcharImgType::New(); m_TissueMask->SetSpacing( m_UpsampledSpacing ); m_TissueMask->SetOrigin( m_Origin ); m_TissueMask->SetDirection( m_DirectionMatrix ); m_TissueMask->SetLargestPossibleRegion( m_UpsampledImageRegion ); m_TissueMask->SetBufferedRegion( m_UpsampledImageRegion ); m_TissueMask->SetRequestedRegion( m_UpsampledImageRegion ); m_TissueMask->Allocate(); m_TissueMask->FillBuffer(1); } // resample frequency map if (m_FrequencyMap.IsNotNull()) { itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetInput(m_FrequencyMap); resampler->SetOutputParametersFromImage(m_FrequencyMap); resampler->SetSize(m_UpsampledImageRegion.GetSize()); resampler->SetOutputSpacing(m_UpsampledSpacing); resampler->Update(); m_FrequencyMap = resampler->GetOutput(); } // initialize volume fraction images m_VolumeFractions.clear(); for (int i=0; iSetSpacing( m_UpsampledSpacing ); tempimg->SetOrigin( m_Origin ); tempimg->SetDirection( m_DirectionMatrix ); tempimg->SetLargestPossibleRegion( m_UpsampledImageRegion ); tempimg->SetBufferedRegion( m_UpsampledImageRegion ); tempimg->SetRequestedRegion( m_UpsampledImageRegion ); tempimg->Allocate(); tempimg->FillBuffer(0); m_VolumeFractions.push_back(tempimg); } // resample fiber bundle for sufficient voxel coverage double segmentVolume = 0.0001; float minSpacing = 1; if(m_UpsampledSpacing[0]GetDeepCopy(); fiberBundle->ResampleFibers(minSpacing/m_VolumeAccuracy); double mmRadius = m_FiberRadius/1000; if (mmRadius>0) segmentVolume = M_PI*mmRadius*mmRadius*minSpacing/m_VolumeAccuracy; // generate double images to work with because we don't want to lose precision // we use a separate image for each compartment model std::vector< DoubleDwiType::Pointer > compartments; for (int i=0; iSetSpacing( m_UpsampledSpacing ); doubleDwi->SetOrigin( m_Origin ); doubleDwi->SetDirection( m_DirectionMatrix ); doubleDwi->SetLargestPossibleRegion( m_UpsampledImageRegion ); doubleDwi->SetBufferedRegion( m_UpsampledImageRegion ); doubleDwi->SetRequestedRegion( m_UpsampledImageRegion ); doubleDwi->SetVectorLength( m_FiberModels[0]->GetNumGradients() ); doubleDwi->Allocate(); DoubleDwiType::PixelType pix; pix.SetSize(m_FiberModels[0]->GetNumGradients()); pix.Fill(0.0); doubleDwi->FillBuffer(pix); compartments.push_back(doubleDwi); } double interpFact = 2*atan(-0.5*m_InterpolationShrink); double maxVolume = 0; MITK_INFO << "Generating signal of " << m_FiberModels.size() << " fiber compartments"; vtkSmartPointer fiberPolyData = fiberBundle->GetFiberPolyData(); boost::progress_display disp(numFibers); for( int i=0; iGetCell(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)); itk::Index<3> idx; itk::ContinuousIndex contIndex; m_TissueMask->TransformPhysicalPointToIndex(vertex, idx); m_TissueMask->TransformPhysicalPointToContinuousIndex(vertex, contIndex); if (!m_UseInterpolation) // use nearest neighbour interpolation { if (!m_TissueMask->GetLargestPossibleRegion().IsInside(idx) || m_TissueMask->GetPixel(idx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(idx); pix += segmentVolume*m_FiberModels[k]->SimulateMeasurement(); doubleDwi->SetPixel(idx, pix ); if (pix[baselineIndex]>maxVolume) maxVolume = pix[baselineIndex]; } 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; } frac_x = atan((0.5-frac_x)*m_InterpolationShrink)/interpFact + 0.5; frac_y = atan((0.5-frac_y)*m_InterpolationShrink)/interpFact + 0.5; frac_z = atan((0.5-frac_z)*m_InterpolationShrink)/interpFact + 0.5; // 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_TissueMask->GetLargestPossibleRegion().IsInside(newIdx) || m_TissueMask->GetPixel(newIdx)<=0) continue; // generate signal for each fiber compartment for (int k=0; kSetFiberDirection(dir); DoubleDwiType::PixelType pix = doubleDwi->GetPixel(newIdx); pix += segmentVolume*frac*m_FiberModels[k]->SimulateMeasurement(); doubleDwi->SetPixel(newIdx, pix ); if (pix[baselineIndex]>maxVolume) maxVolume = pix[baselineIndex]; } } } } } ++disp; } MITK_INFO << "Generating signal of " << m_NonFiberModels.size() << " non-fiber compartments"; ImageRegionIterator it3(m_TissueMask, m_TissueMask->GetLargestPossibleRegion()); boost::progress_display disp3(m_TissueMask->GetLargestPossibleRegion().GetNumberOfPixels()); double voxelVolume = m_UpsampledSpacing[0]*m_UpsampledSpacing[1]*m_UpsampledSpacing[2]; double fact = 1; if (m_FiberRadius<0.0001) fact = voxelVolume/maxVolume; while(!it3.IsAtEnd()) { DoubleDwiType::IndexType index = it3.GetIndex(); if (it3.Get()>0) { // get fiber volume fraction DoubleDwiType::Pointer fiberDwi = compartments.at(0); DoubleDwiType::PixelType fiberPix = fiberDwi->GetPixel(index); // intra axonal compartment if (fact>1) // auto scale intra-axonal if no fiber radius is specified { fiberPix *= fact; fiberDwi->SetPixel(index, fiberPix); } double f = fiberPix[baselineIndex]; if (f>voxelVolume || f>0 && m_EnforcePureFiberVoxels) // more fiber than space in voxel? { fiberDwi->SetPixel(index, fiberPix*voxelVolume/f); - - for (int i=1; iSetPixel(index, pix); - m_VolumeFractions.at(i)->SetPixel(index, 1); - } + m_VolumeFractions.at(0)->SetPixel(index, 1); + +// for (int i=1; iSetPixel(index, pix); +// m_VolumeFractions.at(i)->SetPixel(index, 1); +// } } else { m_VolumeFractions.at(0)->SetPixel(index, f/voxelVolume); double nonf = voxelVolume-f; // non-fiber volume double inter = 0; if (m_FiberModels.size()>1) inter = nonf * f/voxelVolume; // intra-axonal fraction of non fiber compartment scales linearly with f double other = nonf - inter; // rest of compartment double singleinter = inter/(m_FiberModels.size()-1); // adjust non-fiber and intra-axonal signal for (int i=1; iGetPixel(index); if (pix[baselineIndex]>0) pix /= pix[baselineIndex]; pix *= singleinter; doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i)->SetPixel(index, singleinter/voxelVolume); } for (int i=0; iGetPixel(index) + m_NonFiberModels[i]->SimulateMeasurement()*other*m_NonFiberModels[i]->GetWeight(); doubleDwi->SetPixel(index, pix); m_VolumeFractions.at(i+m_FiberModels.size())->SetPixel(index, other/voxelVolume*m_NonFiberModels[i]->GetWeight()); } } } ++it3; ++disp3; } // do k-space stuff DoubleDwiType::Pointer doubleOutImage; if (m_Spikes>0 || m_FrequencyMap.IsNotNull() || !m_KspaceArtifacts.empty() || m_kOffset>0 || m_SimulateRelaxation || m_SimulateEddyCurrents || m_Upsampling>1.00001) { MITK_INFO << "Adjusting complex signal"; doubleOutImage = DoKspaceStuff(compartments); m_SignalScale = 1; } else { MITK_INFO << "Summing compartments"; doubleOutImage = compartments.at(0); for (int i=1; i::Pointer adder = itk::AddImageFilter< DoubleDwiType, DoubleDwiType, DoubleDwiType>::New(); adder->SetInput1(doubleOutImage); adder->SetInput2(compartments.at(i)); adder->Update(); doubleOutImage = adder->GetOutput(); } } MITK_INFO << "Finalizing image"; unsigned int window = 0; unsigned int min = itk::NumericTraits::max(); ImageRegionIterator it4 (outImage, outImage->GetLargestPossibleRegion()); DoubleDwiType::PixelType signal; signal.SetSize(m_FiberModels[0]->GetNumGradients()); boost::progress_display disp4(outImage->GetLargestPossibleRegion().GetNumberOfPixels()); while(!it4.IsAtEnd()) { ++disp4; typename OutputImageType::IndexType index = it4.GetIndex(); signal = doubleOutImage->GetPixel(index)*m_SignalScale; if (m_NoiseModel->GetNoiseVariance() > 0) { DoubleDwiType::PixelType accu = signal; accu.Fill(0.0); for (int i=0; iAddNoise(temp); accu += temp; } signal = accu/m_NumberOfRepetitions; } for (int i=0; i0) signal[i] = floor(signal[i]+0.5); else signal[i] = ceil(signal[i]-0.5); if (!m_FiberModels.at(0)->IsBaselineIndex(i) && signal[i]>window) window = signal[i]; if (!m_FiberModels.at(0)->IsBaselineIndex(i) && signal[i]SetNthOutput(0, outImage); } template< class PixelType > itk::Point TractsToDWIImageFilter< 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 > itk::Vector TractsToDWIImageFilter< 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 TractsToDWIImageFilter< 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 TractsToDWIImageFilter< PixelType >::GetVnlVector(Vector& vector) { vnl_vector_fixed vnlVector; vnlVector[0] = vector[0]; vnlVector[1] = vector[1]; vnlVector[2] = vector[2]; return vnlVector; } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp index e166afc041..749c4f70f8 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp @@ -1,296 +1,302 @@ /*=================================================================== 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 "mitkFiberBundleXReader.h" #include #include #include #include #include #include #include #include #include #include +#include namespace mitk { - void FiberBundleXReader - ::GenerateData() - { +void FiberBundleXReader +::GenerateData() +{ if ( ( ! m_OutputCache ) ) { - Superclass::SetNumberOfRequiredOutputs(0); - this->GenerateOutputInformation(); + Superclass::SetNumberOfRequiredOutputs(0); + this->GenerateOutputInformation(); } if (!m_OutputCache) { - itkWarningMacro("Output cache is empty!"); + itkWarningMacro("Output cache is empty!"); } Superclass::SetNumberOfRequiredOutputs(1); Superclass::SetNthOutput(0, m_OutputCache.GetPointer()); - } +} - void FiberBundleXReader::GenerateOutputInformation() - { +void FiberBundleXReader::GenerateOutputInformation() +{ try { - const std::string& locale = "C"; - const std::string& currLocale = setlocale( LC_ALL, NULL ); - setlocale(LC_ALL, locale.c_str()); - - std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); - ext = itksys::SystemTools::LowerCase(ext); - - vtkSmartPointer chooser=vtkSmartPointer::New(); - chooser->SetFileName(m_FileName.c_str() ); - if( chooser->IsFilePolyData()) - { - MITK_INFO << "Reading vtk fiber bundle"; - vtkSmartPointer reader = vtkSmartPointer::New(); - reader->SetFileName( m_FileName.c_str() ); - reader->Update(); - - if ( reader->GetOutput() != NULL ) - { - vtkSmartPointer fiberPolyData = reader->GetOutput(); + const std::string& locale = "C"; + const std::string& currLocale = setlocale( LC_ALL, NULL ); + setlocale(LC_ALL, locale.c_str()); -// vtkSmartPointer cleaner = vtkSmartPointer::New(); -// cleaner->SetInput(fiberPolyData); -// cleaner->Update(); -// fiberPolyData = cleaner->GetOutput(); + std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); + ext = itksys::SystemTools::LowerCase(ext); - m_OutputCache = OutputType::New(fiberPolyData); + if (ext==".trk") + { + MITK_INFO << "Importing fiber bundle from TrackVis ..."; + m_OutputCache = OutputType::New(); + TrackVis trk; + trk.open(m_FileName); + trk.read(m_OutputCache); + MITK_INFO << "... done"; + return; } - } - else // try to read deprecated fiber bundle file format - { - MITK_INFO << "Reading xml fiber bundle"; - - vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); - vtkSmartPointer cellArray = vtkSmartPointer::New(); - vtkSmartPointer points = vtkSmartPointer::New(); - TiXmlDocument doc( m_FileName ); - if(doc.LoadFile()) + + vtkSmartPointer chooser=vtkSmartPointer::New(); + chooser->SetFileName(m_FileName.c_str() ); + if( chooser->IsFilePolyData()) { - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(0); - - pElem = hDoc.FirstChildElement().Element(); - - // save this for later - hRoot = TiXmlHandle(pElem); - - pElem = hRoot.FirstChildElement("geometry").Element(); - - // read geometry - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - - // read origin - mitk::Point3D origin; - double temp = 0; - pElem->Attribute("origin_x", &temp); - origin[0] = temp; - pElem->Attribute("origin_y", &temp); - origin[1] = temp; - pElem->Attribute("origin_z", &temp); - origin[2] = temp; - geometry->SetOrigin(origin); - - // read spacing - float spacing[3]; - pElem->Attribute("spacing_x", &temp); - spacing[0] = temp; - pElem->Attribute("spacing_y", &temp); - spacing[1] = temp; - pElem->Attribute("spacing_z", &temp); - spacing[2] = temp; - geometry->SetSpacing(spacing); - - // read transform - vtkMatrix4x4* m = vtkMatrix4x4::New(); - pElem->Attribute("xx", &temp); - m->SetElement(0,0,temp); - pElem->Attribute("xy", &temp); - m->SetElement(1,0,temp); - pElem->Attribute("xz", &temp); - m->SetElement(2,0,temp); - pElem->Attribute("yx", &temp); - m->SetElement(0,1,temp); - pElem->Attribute("yy", &temp); - m->SetElement(1,1,temp); - pElem->Attribute("yz", &temp); - m->SetElement(2,1,temp); - pElem->Attribute("zx", &temp); - m->SetElement(0,2,temp); - pElem->Attribute("zy", &temp); - m->SetElement(1,2,temp); - pElem->Attribute("zz", &temp); - m->SetElement(2,2,temp); - - m->SetElement(0,3,origin[0]); - m->SetElement(1,3,origin[1]); - m->SetElement(2,3,origin[2]); - m->SetElement(3,3,1); - geometry->SetIndexToWorldTransformByVtkMatrix(m); - - // read bounds - float bounds[] = {0, 0, 0, 0, 0, 0}; - pElem->Attribute("size_x", &temp); - bounds[1] = temp; - pElem->Attribute("size_y", &temp); - bounds[3] = temp; - pElem->Attribute("size_z", &temp); - bounds[5] = temp; - geometry->SetFloatBounds(bounds); - geometry->SetImageGeometry(true); - - pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); - for( pElem; pElem; pElem=pElem->NextSiblingElement()) - { - TiXmlElement* pElem2 = pElem->FirstChildElement(); - - vtkSmartPointer container = vtkSmartPointer::New(); - - for( pElem2; pElem2; pElem2=pElem2->NextSiblingElement()) - { - itk::Point point; - pElem2->Attribute("pos_x", &temp); - point[0] = temp; - pElem2->Attribute("pos_y", &temp); - point[1] = temp; - pElem2->Attribute("pos_z", &temp); - point[2] = temp; - - geometry->IndexToWorld(point, point); - vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); - container->GetPointIds()->InsertNextId(id); + MITK_INFO << "Reading vtk fiber bundle"; + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName( m_FileName.c_str() ); + reader->Update(); + if ( reader->GetOutput() != NULL ) + { + vtkSmartPointer fiberPolyData = reader->GetOutput(); + m_OutputCache = OutputType::New(fiberPolyData); } - cellArray->InsertNextCell(container); - } - fiberPolyData->SetPoints(points); - fiberPolyData->SetLines(cellArray); - - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInput(fiberPolyData); - cleaner->Update(); - fiberPolyData = cleaner->GetOutput(); - - m_OutputCache = OutputType::New(fiberPolyData); } - else + else // try to read deprecated fiber bundle file format { - MITK_INFO << "could not open xml file"; - throw "could not open xml file"; + MITK_INFO << "Reading xml fiber bundle"; + + vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); + vtkSmartPointer cellArray = vtkSmartPointer::New(); + vtkSmartPointer points = vtkSmartPointer::New(); + TiXmlDocument doc( m_FileName ); + if(doc.LoadFile()) + { + TiXmlHandle hDoc(&doc); + TiXmlElement* pElem; + TiXmlHandle hRoot(0); + + pElem = hDoc.FirstChildElement().Element(); + + // save this for later + hRoot = TiXmlHandle(pElem); + + pElem = hRoot.FirstChildElement("geometry").Element(); + + // read geometry + mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); + + // read origin + mitk::Point3D origin; + double temp = 0; + pElem->Attribute("origin_x", &temp); + origin[0] = temp; + pElem->Attribute("origin_y", &temp); + origin[1] = temp; + pElem->Attribute("origin_z", &temp); + origin[2] = temp; + geometry->SetOrigin(origin); + + // read spacing + float spacing[3]; + pElem->Attribute("spacing_x", &temp); + spacing[0] = temp; + pElem->Attribute("spacing_y", &temp); + spacing[1] = temp; + pElem->Attribute("spacing_z", &temp); + spacing[2] = temp; + geometry->SetSpacing(spacing); + + // read transform + vtkMatrix4x4* m = vtkMatrix4x4::New(); + pElem->Attribute("xx", &temp); + m->SetElement(0,0,temp); + pElem->Attribute("xy", &temp); + m->SetElement(1,0,temp); + pElem->Attribute("xz", &temp); + m->SetElement(2,0,temp); + pElem->Attribute("yx", &temp); + m->SetElement(0,1,temp); + pElem->Attribute("yy", &temp); + m->SetElement(1,1,temp); + pElem->Attribute("yz", &temp); + m->SetElement(2,1,temp); + pElem->Attribute("zx", &temp); + m->SetElement(0,2,temp); + pElem->Attribute("zy", &temp); + m->SetElement(1,2,temp); + pElem->Attribute("zz", &temp); + m->SetElement(2,2,temp); + + m->SetElement(0,3,origin[0]); + m->SetElement(1,3,origin[1]); + m->SetElement(2,3,origin[2]); + m->SetElement(3,3,1); + geometry->SetIndexToWorldTransformByVtkMatrix(m); + + // read bounds + float bounds[] = {0, 0, 0, 0, 0, 0}; + pElem->Attribute("size_x", &temp); + bounds[1] = temp; + pElem->Attribute("size_y", &temp); + bounds[3] = temp; + pElem->Attribute("size_z", &temp); + bounds[5] = temp; + geometry->SetFloatBounds(bounds); + geometry->SetImageGeometry(true); + + pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); + for( pElem; pElem; pElem=pElem->NextSiblingElement()) + { + TiXmlElement* pElem2 = pElem->FirstChildElement(); + + vtkSmartPointer container = vtkSmartPointer::New(); + + for( pElem2; pElem2; pElem2=pElem2->NextSiblingElement()) + { + itk::Point point; + pElem2->Attribute("pos_x", &temp); + point[0] = temp; + pElem2->Attribute("pos_y", &temp); + point[1] = temp; + pElem2->Attribute("pos_z", &temp); + point[2] = temp; + + geometry->IndexToWorld(point, point); + vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); + container->GetPointIds()->InsertNextId(id); + + } + cellArray->InsertNextCell(container); + } + fiberPolyData->SetPoints(points); + fiberPolyData->SetLines(cellArray); + + vtkSmartPointer cleaner = vtkSmartPointer::New(); + cleaner->SetInput(fiberPolyData); + cleaner->Update(); + fiberPolyData = cleaner->GetOutput(); + + m_OutputCache = OutputType::New(fiberPolyData); + } + else + { + MITK_INFO << "could not open xml file"; + throw "could not open xml file"; + } } - } - setlocale(LC_ALL, currLocale.c_str()); - MITK_INFO << "Fiber bundle read"; + setlocale(LC_ALL, currLocale.c_str()); + MITK_INFO << "Fiber bundle read"; } catch(...) { - throw; + throw; } - } +} - void FiberBundleXReader::Update() - { +void FiberBundleXReader::Update() +{ this->GenerateData(); - } +} - const char* FiberBundleXReader - ::GetFileName() const - { +const char* FiberBundleXReader +::GetFileName() const +{ return m_FileName.c_str(); - } +} - void FiberBundleXReader - ::SetFileName(const char* aFileName) - { +void FiberBundleXReader +::SetFileName(const char* aFileName) +{ m_FileName = aFileName; - } +} - const char* FiberBundleXReader - ::GetFilePrefix() const - { +const char* FiberBundleXReader +::GetFilePrefix() const +{ return m_FilePrefix.c_str(); - } +} - void FiberBundleXReader - ::SetFilePrefix(const char* aFilePrefix) - { +void FiberBundleXReader +::SetFilePrefix(const char* aFilePrefix) +{ m_FilePrefix = aFilePrefix; - } +} - const char* FiberBundleXReader - ::GetFilePattern() const - { +const char* FiberBundleXReader +::GetFilePattern() const +{ return m_FilePattern.c_str(); - } +} - void FiberBundleXReader - ::SetFilePattern(const char* aFilePattern) - { +void FiberBundleXReader +::SetFilePattern(const char* aFilePattern) +{ m_FilePattern = aFilePattern; - } +} - bool FiberBundleXReader - ::CanReadFile(const std::string filename, const std::string /*filePrefix*/, const std::string /*filePattern*/) - { +bool FiberBundleXReader +::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 == ".fib") + if (ext == ".fib" || ext == ".trk") { - return true; + return true; } return false; - } +} - BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(const DataObjectIdentifierType &name) - { +BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(const DataObjectIdentifierType &name) +{ itkDebugMacro("MakeOutput(" << name << ")"); if( this->IsIndexedOutputName(name) ) - { - return this->MakeOutput( this->MakeIndexFromOutputName(name) ); - } + { + return this->MakeOutput( this->MakeIndexFromOutputName(name) ); + } return static_cast(OutputType::New().GetPointer()); - } +} - BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(DataObjectPointerArraySizeType /*idx*/) - { +BaseDataSource::DataObjectPointer FiberBundleXReader::MakeOutput(DataObjectPointerArraySizeType /*idx*/) +{ return OutputType::New().GetPointer(); - } +} } //namespace MITK diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp index 93a2a8549f..60b1d0dc4b 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp @@ -1,95 +1,139 @@ /*=================================================================== 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 "mitkFiberBundleXWriter.h" #include #include +#include +#include +#include mitk::FiberBundleXWriter::FiberBundleXWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::FiberBundleXWriter::~FiberBundleXWriter() {} void mitk::FiberBundleXWriter::GenerateData() { try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, NULL ); setlocale(LC_ALL, locale.c_str()); - MITK_INFO << "Writing fiber bundle"; m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to FiberBundleXWriter is NULL!"); return; } else if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } - vtkSmartPointer writer = vtkSmartPointer::New(); - writer->SetInput(input->GetFiberPolyData()); - writer->SetFileName(m_FileName.c_str()); - writer->SetFileTypeToASCII(); - writer->Write(); + + std::string ext = itksys::SystemTools::GetFilenameLastExtension(m_FileName); + + if (ext==".fib" || ext==".vtk") + { + MITK_INFO << "Writing fiber bundle as binary VTK"; + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInput(input->GetFiberPolyData()); + writer->SetFileName(m_FileName.c_str()); + writer->SetFileTypeToBinary(); + writer->Write(); + } + else if (ext==".afib") + { + itksys::SystemTools::ReplaceString(m_FileName,".afib",".fib"); + MITK_INFO << "Writing fiber bundle as ascii VTK"; + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInput(input->GetFiberPolyData()); + writer->SetFileName(m_FileName.c_str()); + writer->SetFileTypeToASCII(); + writer->Write(); + } + else if (ext==".avtk") + { + itksys::SystemTools::ReplaceString(m_FileName,".avtk",".vtk"); + MITK_INFO << "Writing fiber bundle as ascii VTK"; + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInput(input->GetFiberPolyData()); + writer->SetFileName(m_FileName.c_str()); + writer->SetFileTypeToASCII(); + writer->Write(); + } + else if (ext==".trk") + { + MITK_INFO << "Writing fiber bundle as TRK"; + TrackVis trk; + itk::Size<3> size; + size.SetElement(0,input->GetGeometry()->GetExtent(0)); + size.SetElement(1,input->GetGeometry()->GetExtent(1)); + size.SetElement(2,input->GetGeometry()->GetExtent(2)); + trk.create(m_FileName,size,input->GetGeometry()->GetOrigin()); + trk.writeHdr(); + trk.append(input); + } setlocale(LC_ALL, currLocale.c_str()); m_Success = true; MITK_INFO << "Fiber bundle written"; } catch(...) { throw; } } void mitk::FiberBundleXWriter::SetInputFiberBundleX( InputType* diffVolumes ) { this->ProcessObject::SetNthInput( 0, diffVolumes ); } mitk::FiberBundleX* mitk::FiberBundleXWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::FiberBundleXWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".fib"); + possibleFileExtensions.push_back(".afib"); possibleFileExtensions.push_back(".vtk"); + possibleFileExtensions.push_back(".avtk"); + possibleFileExtensions.push_back(".trk"); return possibleFileExtensions; } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h index ae436596a2..3a49873316 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h @@ -1,216 +1,216 @@ /*=================================================================== 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 __mitkFiberBundleXWriter_h #define __mitkFiberBundleXWriter_h #include #include #include "mitkFiberBundleX.h" #include #include namespace mitk { /** * Writes fiber bundles to a file * @ingroup Process */ class FiberTracking_EXPORT FiberBundleXWriter : public mitk::FileWriterWithInformation { public: mitkClassMacro( FiberBundleXWriter, mitk::FileWriterWithInformation ); itkNewMacro( Self ); //mitkWriterMacro; virtual void Write() { if ( this->GetInput() == NULL ) { itkExceptionMacro(<<"Write:Please specify an input!"); return; } /* Fill in image information.*/ this->UpdateOutputInformation(); (*(this->GetInputs().begin()))->SetRequestedRegionToLargestPossibleRegion(); this->PropagateRequestedRegion(NULL); this->UpdateOutputData(NULL); } virtual void Update() { Write(); } typedef mitk::FiberBundleX InputType; /** * Sets the filename of the file to write. * @param FileName the name of the file to write. */ itkSetStringMacro( FileName ); /** * @returns the name of the file to be written to disk. */ itkGetStringMacro( FileName ); /** * @warning multiple write not (yet) supported */ itkSetStringMacro( FilePrefix ); /** * @warning multiple write not (yet) supported */ itkGetStringMacro( FilePrefix ); /** * @warning multiple write not (yet) supported */ itkSetStringMacro( FilePattern ); /** * @warning multiple write not (yet) supported */ itkGetStringMacro( FilePattern ); /** * Sets the input object for the filter. * @param input the diffusion volumes to write to file. */ void SetInputFiberBundleX( InputType* input ); /** * @returns the 0'th input object of the filter. */ InputType* GetInput(); /** * Returns false if an error happened during writing */ itkGetMacro( Success, bool ); /** * @return possible file extensions for the data type associated with the writer */ virtual std::vector GetPossibleFileExtensions(); // FileWriterWithInformation methods virtual const char * GetDefaultFilename() { return "FiberBundle.fib"; } - virtual const char * GetFileDialogPattern() { return "Fiber Bundle (*.fib *.vtk)"; } + virtual const char * GetFileDialogPattern() { return "Fiber Bundle (*.fib *.vtk *.trk *.afib *.avtk)"; } virtual const char * GetDefaultExtension() { return ".fib"; } virtual bool CanWriteBaseDataType(BaseData::Pointer data) { return (dynamic_cast(data.GetPointer()) != NULL); }; virtual void DoWrite(BaseData::Pointer data) { if (CanWriteBaseDataType(data)) { this->SetInputFiberBundleX(dynamic_cast(data.GetPointer())); this->Update(); } }; static const char* XML_GEOMETRY; static const char* XML_MATRIX_XX; static const char* XML_MATRIX_XY; static const char* XML_MATRIX_XZ; static const char* XML_MATRIX_YX; static const char* XML_MATRIX_YY; static const char* XML_MATRIX_YZ; static const char* XML_MATRIX_ZX; static const char* XML_MATRIX_ZY; static const char* XML_MATRIX_ZZ; static const char* XML_ORIGIN_X; static const char* XML_ORIGIN_Y; static const char* XML_ORIGIN_Z; static const char* XML_SPACING_X; static const char* XML_SPACING_Y; static const char* XML_SPACING_Z; static const char* XML_SIZE_X; static const char* XML_SIZE_Y; static const char* XML_SIZE_Z; static const char* XML_FIBER_BUNDLE; static const char* XML_FIBER; static const char* XML_PARTICLE; static const char* XML_ID; static const char* XML_POS_X; static const char* XML_POS_Y; static const char* XML_POS_Z; static const char* VERSION_STRING; static const char* XML_FIBER_BUNDLE_FILE; static const char* XML_FILE_VERSION; static const char* XML_NUM_FIBERS; static const char* XML_NUM_PARTICLES; static const char* ASCII_FILE; static const char* FILE_NAME; protected: FiberBundleXWriter(); virtual ~FiberBundleXWriter(); virtual void GenerateData(); std::string m_FileName; std::string m_FilePrefix; std::string m_FilePattern; bool m_Success; }; } // end of namespace mitk #endif //__mitkFiberBundleXWriter_h diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp new file mode 100644 index 0000000000..77922360e1 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.cpp @@ -0,0 +1,184 @@ +#include + +TrackVis::TrackVis() { filename = ""; fp = NULL; maxSteps = 20000; } + +TrackVis::~TrackVis() { if (fp) fclose( fp ); } + + +// Create a TrackVis file and store standard metadata. The file is ready to append fibers. +// --------------------------------------------------------------------------------------- +short TrackVis::create(string filename , itk::Size<3> size, itk::Point origin) +{ + // prepare the header + for(int i=0; i<3 ;i++) + { + hdr.dim[i] = size.GetElement(i); + hdr.voxel_size[i] = 1; + hdr.origin[i] = 0; + } + hdr.n_scalars = 0; + hdr.n_properties = 0; + sprintf(hdr.voxel_order,"RAS"); + sprintf(hdr.pad2,"LPS"); + hdr.image_orientation_patient[0] = 1.0; + hdr.image_orientation_patient[1] = 0.0; + hdr.image_orientation_patient[2] = 0.0; + hdr.image_orientation_patient[3] = 0.0; + hdr.image_orientation_patient[4] = 1.0; + hdr.image_orientation_patient[5] = 0.0; + hdr.pad1[0] = 0; + hdr.pad1[1] = 0; + hdr.invert_x = 0; + hdr.invert_y = 0; + hdr.invert_z = 0; + hdr.swap_xy = 0; + hdr.swap_yz = 0; + hdr.swap_zx = 0; + hdr.n_count = 0; + hdr.version = 1; + hdr.hdr_size = 1000; + + // write the header to the file + fp = fopen(filename.c_str(),"w+b"); + if (fp == NULL) + { + printf("[ERROR] Unable to create file '%s'\n",filename.c_str()); + return 0; + } + sprintf(hdr.id_string,"TRACK"); + fwrite((char*)&hdr, 1, 1000, fp); + + this->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 TrackVis::open( string filename ) +{ + fp = fopen(filename.c_str(),"r+b"); + if (fp == NULL) + { + printf("[ERROR] Unable to open file '%s'\n",filename.c_str()); + return 0; + } + this->filename = filename; + + return fread((char*)(&hdr), 1, 1000, fp); +} + + + +// Append a fiber to the file +// -------------------------- +short TrackVis::append(mitk::FiberBundleX *fib) +{ + vtkPolyData* poly = fib->GetFiberPolyData(); + for (int i=0; iGetNumFibers(); i++) + { + vtkCell* cell = poly->GetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + unsigned int numSaved, pos = 0; + float tmp[3*maxSteps]; + + if ( numPoints > maxSteps ) + { + printf( "[ERROR] Trying to write a fiber too long!\n" ); + return 0; + } + + numSaved = numPoints; + for(unsigned int i=0; iGetPoint(i); + + tmp[pos++] = -p[0]; + tmp[pos++] = -p[1]; + tmp[pos++] = p[2]; + } + + // write the coordinates to the file + if ( fwrite((char*)&numSaved, 1, 4, fp) != 4 ) + { + printf( "[ERROR] Problems saving the fiber!\n" ); + return 1; + } + if ( fwrite((char*)&tmp, 1, 4*pos, fp) != 4*pos ) + { + printf( "[ERROR] Problems saving the fiber!\n" ); + return 1; + } + } + + return 0; +} + + + +//// Read one fiber from the file +//// ---------------------------- +short TrackVis::read( mitk::FiberBundleX* fib ) +{ + int numPoints; + vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); + vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); + + while (fread((char*)&numPoints, 1, 4, fp)==4) + { + if ( numPoints >= maxSteps || numPoints <= 0 ) + { + printf( "[ERROR] Trying to read a fiber with %d points!\n", numPoints ); + return -1; + } + vtkSmartPointer container = vtkSmartPointer::New(); + + float tmp[3]; + for(int i=0; iInsertNextPoint(tmp); + container->GetPointIds()->InsertNextId(id); + } + vtkNewCells->InsertNextCell(container); + } + + vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); + fiberPolyData->SetPoints(vtkNewPoints); + fiberPolyData->SetLines(vtkNewCells); + fib->SetFiberPolyData(fiberPolyData); + + return numPoints; +} + + + +// Update the field in the header to the new FIBER TOTAL. +// ------------------------------------------------------ +void TrackVis::updateTotal( int totFibers ) +{ + fseek(fp, 1000-12, SEEK_SET); + fwrite((char*)&totFibers, 1, 4, fp); +} + + +void TrackVis::writeHdr() +{ + fseek(fp, 0, SEEK_SET); + fwrite((char*)&hdr, 1, 1000, fp); +} + + +// Close the TrackVis file, but keep the metadata in the header. +// ------------------------------------------------------------- +void TrackVis::close() +{ + fclose(fp); + fp = NULL; +} diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h new file mode 100644 index 0000000000..c7ded97b14 --- /dev/null +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundleX/mitkTrackvis.h @@ -0,0 +1,82 @@ +/*=================================================================== + +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 _TRACKVIS +#define _TRACKVIS + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + +// Structure to hold metadata of a TrackVis file +// --------------------------------------------- +struct TrackVis_header +{ + char id_string[6]; + short int dim[3]; + float voxel_size[3]; + float origin[3]; + short int n_scalars; + char scalar_name[10][20]; + short int n_properties; + char property_name[10][20]; + char reserved[508]; + char voxel_order[4]; + char pad2[4]; + float image_orientation_patient[6]; + char pad1[2]; + unsigned char invert_x; + unsigned char invert_y; + unsigned char invert_z; + unsigned char swap_xy; + unsigned char swap_yz; + unsigned char swap_zx; + int n_count; + int version; + int hdr_size; +}; + +// Class to handle TrackVis files. +// ------------------------------- +class FiberTracking_EXPORT TrackVis +{ +private: + string filename; + FILE* fp; + int maxSteps; // [TODO] should be related to the variable defined for fiber-tracking + +public: + TrackVis_header hdr; + + short create( string filename, itk::Size<3> size, itk::Point origin ); + short open( string filename ); + short read( mitk::FiberBundleX* fib ); + short append( mitk::FiberBundleX* fib ); + void writeHdr(); + void updateTotal( int totFibers ); + void close(); + + TrackVis(); + ~TrackVis(); +}; + +#endif diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberTrackingObjectFactory.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberTrackingObjectFactory.cpp index 442fcbf642..d49ee342c3 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberTrackingObjectFactory.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/mitkFiberTrackingObjectFactory.cpp @@ -1,124 +1,126 @@ #include "mitkFiberTrackingObjectFactory.h" mitk::FiberTrackingObjectFactory::FiberTrackingObjectFactory(bool /*registerSelf*/) :CoreObjectFactoryBase() { static bool alreadyDone = false; if (!alreadyDone) { MITK_DEBUG << "FiberTrackingObjectFactory c'tor" << std::endl; RegisterIOFactories(); mitk::FiberBundleXIOFactory::RegisterOneFactory(); //modernized mitk::FiberBundleXWriterFactory::RegisterOneFactory();//modernized m_FileWriters.push_back( mitk::FiberBundleXWriter::New().GetPointer() );//modernized mitk::CoreObjectFactory::GetInstance()->RegisterExtraFactory(this); CreateFileExtensionsMap(); alreadyDone = true; } } mitk::Mapper::Pointer mitk::FiberTrackingObjectFactory::CreateMapper(mitk::DataNode* node, MapperSlotId id) { mitk::Mapper::Pointer newMapper=NULL; if ( id == mitk::BaseRenderer::Standard2D ) { std::string classname("FiberBundleX"); if(node->GetData() && classname.compare(node->GetData()->GetNameOfClass())==0) { newMapper = mitk::FiberBundleXMapper2D::New(); newMapper->SetDataNode(node); } } else if ( id == mitk::BaseRenderer::Standard3D ) { std::string classname("FiberBundleX"); if(node->GetData() && classname.compare(node->GetData()->GetNameOfClass())==0) { newMapper = mitk::FiberBundleXMapper3D::New(); newMapper->SetDataNode(node); } // classname = "FiberBundleXThreadMonitor"; // if(node->GetData() && classname.compare(node->GetData()->GetNameOfClass())==0) // { // newMapper = mitk::FiberBundleXThreadMonitorMapper3D::New(); // newMapper->SetDataNode(node); // } } return newMapper; } void mitk::FiberTrackingObjectFactory::SetDefaultProperties(mitk::DataNode* node) { std::string classname("FiberBundleX"); if(node->GetData() && classname.compare(node->GetData()->GetNameOfClass())==0) { mitk::FiberBundleXMapper3D::SetDefaultProperties(node); mitk::FiberBundleXMapper2D::SetDefaultProperties(node); } // classname = "FiberBundleXThreadMonitor"; // if(node->GetData() && classname.compare(node->GetData()->GetNameOfClass())==0) // { // mitk::FiberBundleXThreadMonitorMapper3D::SetDefaultProperties(node); // } } const char* mitk::FiberTrackingObjectFactory::GetFileExtensions() { std::string fileExtension; this->CreateFileExtensions(m_FileExtensionsMap, fileExtension); return fileExtension.c_str(); } mitk::CoreObjectFactoryBase::MultimapType mitk::FiberTrackingObjectFactory::GetFileExtensionsMap() { return m_FileExtensionsMap; } const char* mitk::FiberTrackingObjectFactory::GetSaveFileExtensions() { std::string fileExtension; this->CreateFileExtensions(m_SaveFileExtensionsMap, fileExtension); return fileExtension.c_str(); } mitk::CoreObjectFactoryBase::MultimapType mitk::FiberTrackingObjectFactory::GetSaveFileExtensionsMap() { return m_SaveFileExtensionsMap; } void mitk::FiberTrackingObjectFactory::CreateFileExtensionsMap() { m_FileExtensionsMap.insert(std::pair("*.fib", "Fiber Bundle")); m_FileExtensionsMap.insert(std::pair("*.vtk", "Fiber Bundle")); + m_FileExtensionsMap.insert(std::pair("*.trk", "TrackVis Fiber Bundle")); m_SaveFileExtensionsMap.insert(std::pair("*.fib", "Fiber Bundle")); m_SaveFileExtensionsMap.insert(std::pair("*.vtk", "Fiber Bundle")); + m_SaveFileExtensionsMap.insert(std::pair("*.trk", "TrackVis Fiber Bundle")); } void mitk::FiberTrackingObjectFactory::RegisterIOFactories() { } void RegisterFiberTrackingObjectFactory() { static bool oneFiberTrackingObjectFactoryRegistered = false; if ( ! oneFiberTrackingObjectFactoryRegistered ) { MITK_DEBUG << "Registering FiberTrackingObjectFactory..." << std::endl; mitk::CoreObjectFactory::GetInstance()->RegisterExtraFactory(mitk::FiberTrackingObjectFactory::New()); oneFiberTrackingObjectFactoryRegistered = true; } } diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/PeakExtraction.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/PeakExtraction.cpp index c984d03faa..092e723123 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/PeakExtraction.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/PeakExtraction.cpp @@ -1,351 +1,375 @@ /*=================================================================== 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 #include #include #include #include "ctkCommandLineParser.h" #include #include #include #include #include #include - mitk::Image::Pointer LoadData(std::string filename) { if( filename.empty() ) return NULL; const std::string s1="", s2=""; std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( filename, s1, s2, false ); if( infile.empty() ) { MITK_INFO << "File " << filename << " could not be read!"; return NULL; } mitk::BaseData::Pointer baseData = infile.at(0); return dynamic_cast(baseData.GetPointer()); } -int PeakExtraction(int argc, char* argv[]) + +template +int Start(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("image", "i", ctkCommandLineParser::String, "sh coefficient image", us::Any(), false); parser.addArgument("outroot", "o", ctkCommandLineParser::String, "output root", us::Any(), false); parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "0=no norm, 1=max norm, 2=single vec norm", 1, true); parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "maximum number of extracted peaks", 2, true); parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "peak threshold relative to largest peak", 0.4, true); parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "absolute peak threshold weighted with local GFA value", 0.06, true); parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "do not flip input image to match MITK coordinate convention"); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; // mandatory arguments string imageName = us::any_cast(parsedArgs["image"]); string outRoot = us::any_cast(parsedArgs["outroot"]); // optional arguments string maskImageName(""); if (parsedArgs.count("mask")) maskImageName = us::any_cast(parsedArgs["mask"]); int normalization = 1; if (parsedArgs.count("normalization")) normalization = us::any_cast(parsedArgs["normalization"]); int numPeaks = 2; if (parsedArgs.count("numpeaks")) numPeaks = us::any_cast(parsedArgs["numpeaks"]); float peakThres = 0.4; if (parsedArgs.count("peakthres")) peakThres = us::any_cast(parsedArgs["peakthres"]); float absPeakThres = 0.06; if (parsedArgs.count("abspeakthres")) absPeakThres = us::any_cast(parsedArgs["abspeakthres"]); bool noFlip = false; if (parsedArgs.count("noFlip")) noFlip = us::any_cast(parsedArgs["noFlip"]); MITK_INFO << "image: " << imageName; MITK_INFO << "outroot: " << outRoot; if (!maskImageName.empty()) MITK_INFO << "mask: " << maskImageName; else MITK_INFO << "no mask image selected"; MITK_INFO << "numpeaks: " << numPeaks; MITK_INFO << "peakthres: " << peakThres; MITK_INFO << "abspeakthres: " << absPeakThres; + MITK_INFO << "shOrder: " << shOrder; try { mitk::CoreObjectFactory::GetInstance(); RegisterDiffusionCoreObjectFactory(); mitk::Image::Pointer image = LoadData(imageName); mitk::Image::Pointer mask = LoadData(maskImageName); typedef itk::Image ItkUcharImgType; - typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, 4, 20242 > MaximaExtractionFilterType; - MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); + typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; + typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); int toolkitConvention = 0; if (parsedArgs.count("shConvention")) { string convention = us::any_cast(parsedArgs["shConvention"]).c_str(); if ( boost::algorithm::equals(convention, "FSL") ) { toolkitConvention = 1; MITK_INFO << "Using FSL SH-basis"; } else if ( boost::algorithm::equals(convention, "MRtrix") ) { toolkitConvention = 2; MITK_INFO << "Using MRtrix SH-basis"; } else MITK_INFO << "Using MITK SH-basis"; } else MITK_INFO << "Using MITK SH-basis"; ItkUcharImgType::Pointer itkMaskImage = NULL; if (mask.IsNotNull()) { try{ itkMaskImage = ItkUcharImgType::New(); mitk::CastToItkImage(mask, itkMaskImage); filter->SetMaskImage(itkMaskImage); } catch(...) { } } if (toolkitConvention>0) { MITK_INFO << "Converting coefficient image to MITK format"; - typedef itk::ShCoefficientImageImporter< float, 4 > ConverterType; + typedef itk::ShCoefficientImageImporter< float, shOrder > ConverterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); itk::Image< float, 4 >::Pointer itkImage = caster->GetOutput(); - ConverterType::Pointer converter = ConverterType::New(); + typename ConverterType::Pointer converter = ConverterType::New(); if (noFlip) { converter->SetInputImage(itkImage); } else { + MITK_INFO << "Flipping image"; itk::FixedArray flipAxes; flipAxes[0] = true; flipAxes[1] = true; flipAxes[2] = false; flipAxes[3] = false; itk::FlipImageFilter< itk::Image< float, 4 > >::Pointer flipper = itk::FlipImageFilter< itk::Image< float, 4 > >::New(); flipper->SetInput(itkImage); flipper->SetFlipAxes(flipAxes); flipper->Update(); itk::Image< float, 4 >::Pointer flipped = flipper->GetOutput(); flipped->SetDirection(itkImage->GetDirection()); flipped->SetOrigin(itkImage->GetOrigin()); converter->SetInputImage(flipped); } + MITK_INFO << "Starting conversion"; switch (toolkitConvention) { case 1: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; case 2: converter->SetToolkit(ConverterType::MRTRIX); filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); break; default: converter->SetToolkit(ConverterType::FSL); filter->SetToolkit(MaximaExtractionFilterType::FSL); break; } converter->GenerateData(); filter->SetInput(converter->GetCoefficientImage()); - - // write qbi - ConverterType::QballImageType::Pointer itkQbi = converter->GetQballImage(); - - if (itkMaskImage.IsNotNull()) - { - itkQbi->SetDirection(itkMaskImage->GetDirection()); - itkQbi->SetOrigin(itkMaskImage->GetOrigin()); - } - - mitk::QBallImage::Pointer mitkQbi = mitk::QBallImage::New(); - mitkQbi->InitializeByItk( itkQbi.GetPointer() ); - mitkQbi->SetVolume( itkQbi->GetBufferPointer() ); - - string outfilename = outRoot; - outfilename.append("_QBI.qbi"); - MITK_INFO << "writing " << outfilename; - - mitk::NrrdQBallImageWriter::Pointer writer = mitk::NrrdQBallImageWriter::New(); - writer->SetFileName(outfilename.c_str()); - writer->DoWrite(mitkQbi.GetPointer()); } else { try{ - typedef mitk::ImageToItk< MaximaExtractionFilterType::CoefficientImageType > CasterType; - CasterType::Pointer caster = CasterType::New(); + typedef mitk::ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; + typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); } catch(...) { MITK_INFO << "wrong image type"; return EXIT_FAILURE; } } filter->SetMaxNumPeaks(numPeaks); filter->SetPeakThreshold(peakThres); filter->SetAbsolutePeakThreshold(absPeakThres); filter->SetAngularThreshold(1); switch (normalization) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } + MITK_INFO << "Starting extraction"; filter->Update(); // write direction images { - typedef MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; - ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); + typedef typename MaximaExtractionFilterType::ItkDirectionImageContainer ItkDirectionImageContainer; + typename ItkDirectionImageContainer::Pointer container = filter->GetDirectionImageContainer(); for (int i=0; iSize(); i++) { - MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); + typename MaximaExtractionFilterType::ItkDirectionImage::Pointer itkImg = container->GetElement(i); if (itkMaskImage.IsNotNull()) { itkImg->SetDirection(itkMaskImage->GetDirection()); itkImg->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot; outfilename.append("_DIRECTION_"); outfilename.append(boost::lexical_cast(i)); outfilename.append(".nrrd"); MITK_INFO << "writing " << outfilename; - typedef itk::ImageFileWriter< MaximaExtractionFilterType::ItkDirectionImage > WriterType; - WriterType::Pointer writer = WriterType::New(); + typedef itk::ImageFileWriter< typename MaximaExtractionFilterType::ItkDirectionImage > WriterType; + typename WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(itkImg); writer->Update(); } } // write num directions image { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); if (itkMaskImage.IsNotNull()) { numDirImage->SetDirection(itkMaskImage->GetDirection()); numDirImage->SetOrigin(itkMaskImage->GetOrigin()); } string outfilename = outRoot.c_str(); outfilename.append("_NUM_DIRECTIONS.nrrd"); MITK_INFO << "writing " << outfilename; typedef itk::ImageFileWriter< ItkUcharImgType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outfilename); writer->SetInput(numDirImage); writer->Update(); } // write vector field { mitk::FiberBundleX::Pointer directions = filter->GetOutputFiberBundle(); string outfilename = outRoot.c_str(); outfilename.append("_VECTOR_FIELD.fib"); mitk::FiberBundleXWriter::Pointer fibWriter = mitk::FiberBundleXWriter::New(); fibWriter->SetFileName(outfilename.c_str()); fibWriter->DoWrite(directions.GetPointer()); } } 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; } + +int PeakExtraction(int argc, char* argv[]) +{ + ctkCommandLineParser parser; + parser.setArgumentPrefix("--", "-"); + parser.addArgument("image", "i", ctkCommandLineParser::String, "sh coefficient image", us::Any(), false); + parser.addArgument("shOrder", "sh", ctkCommandLineParser::Int, "spherical harmonics order"); + parser.addArgument("outroot", "o", ctkCommandLineParser::String, "output root", us::Any(), false); + parser.addArgument("mask", "m", ctkCommandLineParser::String, "mask image"); + parser.addArgument("normalization", "n", ctkCommandLineParser::Int, "0=no norm, 1=max norm, 2=single vec norm", 1, true); + parser.addArgument("numpeaks", "p", ctkCommandLineParser::Int, "maximum number of extracted peaks", 2, true); + parser.addArgument("peakthres", "r", ctkCommandLineParser::Float, "peak threshold relative to largest peak", 0.4, true); + parser.addArgument("abspeakthres", "a", ctkCommandLineParser::Float, "absolute peak threshold weighted with local GFA value", 0.06, true); + parser.addArgument("shConvention", "s", ctkCommandLineParser::String, "use specified SH-basis (MITK, FSL, MRtrix)", string("MITK"), true); + parser.addArgument("noFlip", "f", ctkCommandLineParser::Bool, "do not flip input image to match MITK coordinate convention"); + + map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; + + + int shOrder = -1; + if (parsedArgs.count("shOrder")) + shOrder = us::any_cast(parsedArgs["shOrder"]); + + switch (shOrder) + { + case 4: + return Start<4>(argc, argv); + case 6: + return Start<6>(argc, argv); + case 8: + return Start<8>(argc, argv); + case 10: + return Start<10>(argc, argv); + case 12: + return Start<12>(argc, argv); + } + return EXIT_FAILURE; +} RegisterFiberTrackingMiniApp(PeakExtraction); diff --git a/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp index a704a91c2d..adb943a9c7 100755 --- a/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp +++ b/Modules/DiffusionImaging/FiberTracking/MiniApps/StreamlineTracking.cpp @@ -1,174 +1,174 @@ /*=================================================================== 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 #include "ctkCommandLineParser.h" int StreamlineTracking(int argc, char* argv[]) { ctkCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", ctkCommandLineParser::String, "input tensor image (.dti)", us::Any(), false); - parser.addArgument("seed", "s", ctkCommandLineParser::String, "binary seed image"); - parser.addArgument("mask", "m", ctkCommandLineParser::String, "binary mask image"); - parser.addArgument("minFA", "t", ctkCommandLineParser::Float, "minimum fractional anisotropy threshold", 0.15, true); + parser.addArgument("seed", "si", ctkCommandLineParser::String, "binary seed image"); + parser.addArgument("mask", "mi", ctkCommandLineParser::String, "binary mask image"); + parser.addArgument("minFA", "fa", ctkCommandLineParser::Float, "minimum fractional anisotropy threshold", 0.15, true); parser.addArgument("minCurv", "c", ctkCommandLineParser::Float, "minimum curvature radius in mm (default = 0.5*minimum-spacing)"); parser.addArgument("stepSize", "s", ctkCommandLineParser::Float, "stepsize in mm (default = 0.1*minimum-spacing)"); parser.addArgument("tendf", "f", ctkCommandLineParser::Float, "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0, true); parser.addArgument("tendg", "g", ctkCommandLineParser::Float, "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0, true); parser.addArgument("numSeeds", "n", ctkCommandLineParser::Int, "Number of seeds per voxel.", 1, true); parser.addArgument("minLength", "l", ctkCommandLineParser::Float, "minimum fiber length in mm", 20, true); - parser.addArgument("interpolate", "a", ctkCommandLineParser::Bool, "Use linear interpolation", false, true); + parser.addArgument("interpolate", "ip", ctkCommandLineParser::Bool, "Use linear interpolation", false, true); parser.addArgument("outFile", "o", ctkCommandLineParser::String, "output fiber bundle (.fib)", us::Any(), false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string dtiFileName = us::any_cast(parsedArgs["input"]); string outFileName = us::any_cast(parsedArgs["outFile"]); float minFA = 0.15; float minCurv = -1; float stepSize = -1; float tendf = 1; float tendg = 0; float minLength = 20; int numSeeds = 1; bool interpolate = false; if (parsedArgs.count("minCurv")) minCurv = us::any_cast(parsedArgs["minCurv"]); if (parsedArgs.count("minFA")) minFA = us::any_cast(parsedArgs["minFA"]); if (parsedArgs.count("stepSize")) stepSize = us::any_cast(parsedArgs["stepSize"]); if (parsedArgs.count("tendf")) tendf = us::any_cast(parsedArgs["tendf"]); if (parsedArgs.count("tendg")) tendg = us::any_cast(parsedArgs["tendg"]); if (parsedArgs.count("minLength")) minLength = us::any_cast(parsedArgs["minLength"]); if (parsedArgs.count("numSeeds")) numSeeds = us::any_cast(parsedArgs["numSeeds"]); if (parsedArgs.count("interpolate")) interpolate = us::any_cast(parsedArgs["interpolate"]); try { RegisterDiffusionCoreObjectFactory(); RegisterFiberTrackingObjectFactory(); // load input image const std::string s1="", s2=""; std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( dtiFileName, s1, s2, false ); MITK_INFO << "Loading tensor image ..."; typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(infile.at(0).GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); MITK_INFO << "Loading seed image ..."; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; mitk::Image::Pointer mitkSeedImage = NULL; if (parsedArgs.count("seed")) mitkSeedImage = mitk::IOUtil::LoadImage(us::any_cast(parsedArgs["seed"])); MITK_INFO << "Loading mask image ..."; mitk::Image::Pointer mitkMaskImage = NULL; if (parsedArgs.count("mask")) mitkMaskImage = mitk::IOUtil::LoadImage(us::any_cast(parsedArgs["mask"])); // instantiate tracker typedef itk::StreamlineTrackingFilter< float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(itk_dti); filter->SetSeedsPerVoxel(numSeeds); filter->SetFaThreshold(minFA); filter->SetMinCurvatureRadius(minCurv); filter->SetStepSize(stepSize); filter->SetF(tendf); filter->SetG(tendg); filter->SetInterpolate(interpolate); filter->SetMinTractLength(minLength); if (mitkSeedImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkSeedImage, mask); filter->SetSeedImage(mask); } if (mitkMaskImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMaskImage, mask); filter->SetMaskImage(mask); } filter->Update(); vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); if ( fiberBundle->GetNumberOfLines()==0 ) { MITK_INFO << "No fibers reconstructed. Check parametrization."; return EXIT_FAILURE; } mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); mitk::CoreObjectFactory::FileWriterList fileWriters = mitk::CoreObjectFactory::GetInstance()->GetFileWriters(); for (mitk::CoreObjectFactory::FileWriterList::iterator it = fileWriters.begin() ; it != fileWriters.end() ; ++it) { if ( (*it)->CanWriteBaseDataType(fib.GetPointer()) ) { (*it)->SetFileName( outFileName.c_str() ); (*it)->DoWrite( fib.GetPointer() ); } } } 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(StreamlineTracking); diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index ab9a5b3c39..d29b958166 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,106 +1,108 @@ 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/mitkTrackvis.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/FiberBundleX/mitkTrackvis.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 ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppImageProcessingPerspective.cpp b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppImageProcessingPerspective.cpp index 4a4a4bfc94..3c187497b1 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppImageProcessingPerspective.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppImageProcessingPerspective.cpp @@ -1,56 +1,56 @@ /*=================================================================== 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 "QmitkDIAppImageProcessingPerspective.h" #include "berryIViewLayout.h" void QmitkDIAppImageProcessingPerspective::CreateInitialLayout(berry::IPageLayout::Pointer layout) { ///////////////////////////////////////////////////// // all di-app perspectives should have the following: ///////////////////////////////////////////////////// std::string editorArea = layout->GetEditorArea(); layout->AddStandaloneView("org.mitk.views.datamanager", false, berry::IPageLayout::LEFT, 0.3f, editorArea); layout->AddStandaloneView("org.mitk.views.controlvisualizationpropertiesview", false, berry::IPageLayout::BOTTOM, 0.15f, "org.mitk.views.datamanager"); berry::IFolderLayout::Pointer left = layout->CreateFolder("org.mbi.diffusionimaginginternal.leftcontrols", berry::IPageLayout::BOTTOM, 0.1f, "org.mitk.views.controlvisualizationpropertiesview"); layout->AddStandaloneViewPlaceholder("org.mitk.views.imagenavigator", berry::IPageLayout::BOTTOM, .4f, "org.mbi.diffusionimaginginternal.leftcontrols", false); ///////////////////////////////////////////// // here goes the perspective specific stuff ///////////////////////////////////////////// left->AddView("org.mitk.views.segmentation"); berry::IViewLayout::Pointer lo = layout->GetViewLayout("org.mitk.views.segmentation"); lo->SetCloseable(false); left->AddView("org.mitk.views.segmentationboolean"); - lo = layout->GetViewLayout("org.mitk.views.segmentationutilities"); + lo = layout->GetViewLayout("org.mitk.views.segmentationboolean"); lo->SetCloseable(false); left->AddView("org.mitk.views.basicimageprocessing"); lo = layout->GetViewLayout("org.mitk.views.basicimageprocessing"); lo->SetCloseable(false); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppUtilityPerspective.cpp b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppUtilityPerspective.cpp index cf3e502334..b019aee55e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppUtilityPerspective.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/internal/Perspectives/QmitkDIAppUtilityPerspective.cpp @@ -1,52 +1,52 @@ /*=================================================================== 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 "QmitkDIAppUtilityPerspective.h" #include "berryIViewLayout.h" void QmitkDIAppUtilityPerspective::CreateInitialLayout(berry::IPageLayout::Pointer layout) { ///////////////////////////////////////////////////// // all di-app perspectives should have the following: ///////////////////////////////////////////////////// std::string editorArea = layout->GetEditorArea(); layout->AddStandaloneView("org.mitk.views.datamanager", false, berry::IPageLayout::LEFT, 0.3f, editorArea); layout->AddStandaloneView("org.mitk.views.controlvisualizationpropertiesview", false, berry::IPageLayout::BOTTOM, .15f, "org.mitk.views.datamanager"); berry::IFolderLayout::Pointer left = layout->CreateFolder("org.mbi.diffusionimaginginternal.leftcontrols", berry::IPageLayout::BOTTOM, 0.1f, "org.mitk.views.controlvisualizationpropertiesview"); layout->AddStandaloneViewPlaceholder("org.mitk.views.imagenavigator", berry::IPageLayout::BOTTOM, .4f, "org.mbi.diffusionimaginginternal.leftcontrols", false); ///////////////////////////////////////////// // here goes the perspective specific stuff ///////////////////////////////////////////// - left->AddView("org.mitk.views.propertylistview"); + left->AddView("org.mitk.views.properties"); berry::IViewLayout::Pointer lo = layout->GetViewLayout("org.mitk.views.properties"); lo->SetCloseable(false); left->AddView("org.blueberry.views.logview"); lo = layout->GetViewLayout("org.blueberry.views.logview"); lo->SetCloseable(false); }