Skip to content

Commit

Permalink
Add tensor tool for flipping tensor axes.
Browse files Browse the repository at this point in the history
  • Loading branch information
astamm committed May 31, 2018
1 parent b226b68 commit f7e9e03
Show file tree
Hide file tree
Showing 5 changed files with 256 additions and 0 deletions.
1 change: 1 addition & 0 deletions Anima/diffusion/dti_tools/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
add_subdirectory(compute_dti_scalar_maps)
add_subdirectory(flip_tensors)
36 changes: 36 additions & 0 deletions Anima/diffusion/dti_tools/flip_tensors/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -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()
63 changes: 63 additions & 0 deletions Anima/diffusion/dti_tools/flip_tensors/animaFlipTensors.cxx
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
#include <animaFlipTensorsImageFilter.h>
#include <tclap/CmdLine.h>
#include <animaReadWriteFunctions.h>

//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)<<"%"<<std::flush;
}

int main(int argc, char **argv)
{

TCLAP::CmdLine cmd("Flip tensors in a DTI volume.\nINRIA / IRISA - VisAGeS Team",' ',ANIMA_VERSION);

TCLAP::ValueArg<std::string> inArg("i", "input", "Input tensor image", true, "", "input tensor image", cmd);
TCLAP::ValueArg<std::string> outArg("o", "output", "Output flipped tensor image", true, "", "output image", cmd);

TCLAP::ValueArg<std::string> 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<unsigned int> 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<InputImageType>(inArg.getValue()));
filter->SetFlipXAxis(xDirArg.isSet());
filter->SetFlipYAxis(yDirArg.isSet());
filter->SetFlipZAxis(zDirArg.isSet());

if (maskArg.getValue() != "")
filter->SetComputationMask(anima::readImage<MaskImageType>(maskArg.getValue()));

filter->SetNumberOfThreads(nbpArg.getValue());
filter->AddObserver(itk::ProgressEvent(), callback );
filter->Update();
std::cout << std::endl;

anima::writeImage<OutputImageType>(outArg.getValue(), filter->GetOutput());

return EXIT_SUCCESS;
}// end of main
Original file line number Diff line number Diff line change
@@ -0,0 +1,70 @@
#pragma once

#include <animaMaskedImageToImageFilter.h>
#include <itkVectorImage.h>

namespace anima
{

template <unsigned int ImageDimension = 3>
class FlipTensorsImageFilter :
public anima::MaskedImageToImageFilter< itk::VectorImage <float, ImageDimension>,
itk::VectorImage <float, ImageDimension> >
{
public:
/** Standard class typedefs. */
typedef FlipTensorsImageFilter Self;

typedef itk::VectorImage<float,ImageDimension> InputImageType;
typedef typename InputImageType::PixelType InputPixelType;

typedef itk::VectorImage<float,ImageDimension> OutputImageType;
typedef typename OutputImageType::PixelType OutputPixelType;

typedef anima::MaskedImageToImageFilter<InputImageType, OutputImageType> Superclass;

typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> 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"
Original file line number Diff line number Diff line change
@@ -0,0 +1,86 @@
#pragma once

#include <animaBaseTensorTools.h>
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionIterator.h>

namespace anima
{

template <unsigned int ImageDimension>
void
FlipTensorsImageFilter<ImageDimension>
::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 <unsigned int ImageDimension>
void
FlipTensorsImageFilter<ImageDimension>
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
itk::ThreadIdType threadId)
{
itk::ImageRegionConstIterator<InputImageType> inItr(this->GetInput(), outputRegionForThread);
itk::ImageRegionConstIterator<MaskImageType> maskItr(this->GetComputationMask(), outputRegionForThread);
itk::ImageRegionIterator<OutputImageType> outItr(this->GetOutput(), outputRegionForThread);

typedef vnl_matrix<double> MatrixType;
typedef vnl_diag_matrix<double> EigenVectorType;

EigenVectorType eVals(3);
MatrixType eVecs(3,3);
MatrixType tensorSymMatrix(3,3);
InputPixelType inTensor(6);
OutputPixelType outTensor(6);

itk::SymmetricEigenAnalysis<MatrixType, EigenVectorType, MatrixType> 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

0 comments on commit f7e9e03

Please sign in to comment.