diff --git a/Anima/diffusion/dti_tools/CMakeLists.txt b/Anima/diffusion/dti_tools/CMakeLists.txt index 8ab1534e9..7c9e02c31 100644 --- a/Anima/diffusion/dti_tools/CMakeLists.txt +++ b/Anima/diffusion/dti_tools/CMakeLists.txt @@ -1 +1,2 @@ add_subdirectory(compute_dti_scalar_maps) +add_subdirectory(flip_tensors) diff --git a/Anima/diffusion/dti_tools/flip_tensors/CMakeLists.txt b/Anima/diffusion/dti_tools/flip_tensors/CMakeLists.txt new file mode 100644 index 000000000..017f84707 --- /dev/null +++ b/Anima/diffusion/dti_tools/flip_tensors/CMakeLists.txt @@ -0,0 +1,36 @@ +if(BUILD_TOOLS) + +project(animaFlipTensors) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + +## ############################################################################# +## add executable +## ############################################################################# + +add_executable(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + ${ITKIO_LIBRARIES} + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensors.cxx b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensors.cxx new file mode 100644 index 000000000..72d2e8580 --- /dev/null +++ b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensors.cxx @@ -0,0 +1,63 @@ +#include +#include +#include + +//Update progression of the process +void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData) +{ + itk::ProcessObject * processObject = (itk::ProcessObject*) caller; + std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"< inArg("i", "input", "Input tensor image", true, "", "input tensor image", cmd); + TCLAP::ValueArg outArg("o", "output", "Output flipped tensor image", true, "", "output image", cmd); + + TCLAP::ValueArg maskArg("m", "mask", "Computation mask", false, "", "mask", cmd); + + TCLAP::SwitchArg xDirArg("X", "xdir", "Flip X axis", cmd, false); + TCLAP::SwitchArg yDirArg("Y", "ydir", "Flip Y axis", cmd, false); + TCLAP::SwitchArg zDirArg("Z", "zdir", "Flip Z axis", cmd, false); + + TCLAP::ValueArg nbpArg("p", "nthreads", "Number of thread to use (default: all)", false, itk::MultiThreader::GetGlobalDefaultNumberOfThreads(), "number of thread", cmd); + + try + { + cmd.parse(argc,argv); + } + catch (TCLAP::ArgException& e) + { + std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; + return EXIT_FAILURE; + } + + itk::CStyleCommand::Pointer callback = itk::CStyleCommand::New(); + callback->SetCallback(eventCallback); + + typedef anima::FlipTensorsImageFilter<3> FilterType; + typedef FilterType::InputImageType InputImageType; + typedef FilterType::OutputImageType OutputImageType; + typedef FilterType::MaskImageType MaskImageType; + + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(anima::readImage(inArg.getValue())); + filter->SetFlipXAxis(xDirArg.isSet()); + filter->SetFlipYAxis(yDirArg.isSet()); + filter->SetFlipZAxis(zDirArg.isSet()); + + if (maskArg.getValue() != "") + filter->SetComputationMask(anima::readImage(maskArg.getValue())); + + filter->SetNumberOfThreads(nbpArg.getValue()); + filter->AddObserver(itk::ProgressEvent(), callback ); + filter->Update(); + std::cout << std::endl; + + anima::writeImage(outArg.getValue(), filter->GetOutput()); + + return EXIT_SUCCESS; +}// end of main diff --git a/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.h b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.h new file mode 100644 index 000000000..d17e3c3a2 --- /dev/null +++ b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.h @@ -0,0 +1,70 @@ +#pragma once + +#include +#include + +namespace anima +{ + +template +class FlipTensorsImageFilter : + public anima::MaskedImageToImageFilter< itk::VectorImage , + itk::VectorImage > +{ +public: + /** Standard class typedefs. */ + typedef FlipTensorsImageFilter Self; + + typedef itk::VectorImage InputImageType; + typedef typename InputImageType::PixelType InputPixelType; + + typedef itk::VectorImage OutputImageType; + typedef typename OutputImageType::PixelType OutputPixelType; + + typedef anima::MaskedImageToImageFilter Superclass; + + typedef itk::SmartPointer Pointer; + typedef itk::SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Run-time type information (and related methods) */ + itkTypeMacro(FlipTensorsImageFilter, MaskedImageToImageFilter) + + /** Superclass typedefs. */ + typedef typename Superclass::InputImageRegionType InputImageRegionType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + typedef typename Superclass::MaskImageType MaskImageType; + + typedef typename InputImageType::Pointer InputPointerType; + typedef typename OutputImageType::Pointer OutputPointerType; + + itkSetMacro(FlipXAxis, bool) + itkSetMacro(FlipYAxis, bool) + itkSetMacro(FlipZAxis, bool) + +protected: + FlipTensorsImageFilter() : Superclass() + { + m_FlipXAxis = false; + m_FlipYAxis = false; + m_FlipZAxis = false; + } + virtual ~FlipTensorsImageFilter() {} + + void GenerateOutputInformation() ITK_OVERRIDE; + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + itk::ThreadIdType threadId) ITK_OVERRIDE; + +private: + ITK_DISALLOW_COPY_AND_ASSIGN(FlipTensorsImageFilter); + + bool m_FlipXAxis, m_FlipYAxis, m_FlipZAxis; + + static const unsigned int m_NumberOfComponents = 6; +}; + +} // end of namespace anima + +#include "animaFlipTensorsImageFilter.hxx" diff --git a/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.hxx b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.hxx new file mode 100644 index 000000000..9d68350ce --- /dev/null +++ b/Anima/diffusion/dti_tools/flip_tensors/animaFlipTensorsImageFilter.hxx @@ -0,0 +1,86 @@ +#pragma once + +#include +#include +#include + +namespace anima +{ + +template +void +FlipTensorsImageFilter +::GenerateOutputInformation() +{ + // Override the method in itkImageSource, so we can set the vector length of + // the output itk::VectorImage + + this->Superclass::GenerateOutputInformation(); + + OutputImageType *output = this->GetOutput(); + output->SetVectorLength(m_NumberOfComponents); +} + +template +void +FlipTensorsImageFilter +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + itk::ThreadIdType threadId) +{ + itk::ImageRegionConstIterator inItr(this->GetInput(), outputRegionForThread); + itk::ImageRegionConstIterator maskItr(this->GetComputationMask(), outputRegionForThread); + itk::ImageRegionIterator outItr(this->GetOutput(), outputRegionForThread); + + typedef vnl_matrix MatrixType; + typedef vnl_diag_matrix EigenVectorType; + + EigenVectorType eVals(3); + MatrixType eVecs(3,3); + MatrixType tensorSymMatrix(3,3); + InputPixelType inTensor(6); + OutputPixelType outTensor(6); + + itk::SymmetricEigenAnalysis eigenAnalysis(3); + + while (!maskItr.IsAtEnd()) + { + if (maskItr.Get() == 0) + { + outTensor.Fill(0.0); + outItr.Set(outTensor); + + ++inItr; + ++outItr; + ++maskItr; + continue; + } + + inTensor = inItr.Get(); + anima::GetTensorFromVectorRepresentation(inTensor, tensorSymMatrix,3,false); + + eigenAnalysis.ComputeEigenValuesAndVectors(tensorSymMatrix, eVals, eVecs); + + if (m_FlipXAxis) + for (unsigned int i = 0;i < 3;++i) + eVecs(i,0) *= -1.0; + + if (m_FlipYAxis) + for (unsigned int i = 0;i < 3;++i) + eVecs(i,1) *= -1.0; + + if (m_FlipZAxis) + for (unsigned int i = 0;i < 3;++i) + eVecs(i,2) *= -1.0; + + anima::RecomposeTensor(eVals, eVecs, tensorSymMatrix); + anima::GetVectorRepresentation(tensorSymMatrix,outTensor,m_NumberOfComponents,false); + + outItr.Set(outTensor); + + ++inItr; + ++outItr; + ++maskItr; + } +} + +} // end of namespace anima