diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp index b0212be840..2dbeafe76a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkBatchedRegistration.cpp @@ -1,117 +1,155 @@ #include "mitkBatchedRegistration.h" #include "mitkPyramidImageRegistrationMethod.h" #include "mitkDiffusionImage.h" #include +#include +#include "itkB0ImageExtractionImageFilter.h" #include #include // VTK #include mitk::BatchedRegistration::BatchedRegistration() : m_RegisteredImagesValid(false) { } void mitk::BatchedRegistration::SetFixedImage(mitk::Image::Pointer& fixedImage) { m_FixedImage = fixedImage; } void mitk::BatchedRegistration::SetMovingReferenceImage(Image::Pointer &movingImage) { m_MovingReference = movingImage; m_RegisteredImagesValid = false; } void mitk::BatchedRegistration::SetBatch(std::vector imageBatch) { m_ImageBatch.clear(); m_ImageBatch = imageBatch; } std::vector mitk::BatchedRegistration::GetRegisteredImages() { if (!m_RegisteredImagesValid) { m_RegisteredImages.clear(); // First transform moving reference image TransformType transf = GetTransformation(m_FixedImage, m_MovingReference); // store it as first element in vector ApplyTransformationToImage(m_MovingReference,transf); m_RegisteredImages.push_back(m_MovingReference); // apply transformation to whole batch std::vector::const_iterator itEnd = m_ImageBatch.end(); for (std::vector::iterator it = m_ImageBatch.begin(); it != itEnd; ++it) { ApplyTransformationToImage(*it,transf); m_RegisteredImages.push_back(*it); } } return m_RegisteredImages; } void mitk::BatchedRegistration::ApplyTransformationToImage(mitk::Image::Pointer &img, const mitk::BatchedRegistration::TransformType &transformation) const { itk::ScalableAffineTransform::Pointer rotationTransform; vtkMatrix4x4* transformationMatrix = vtkMatrix4x4::New(); + vnl_matrix_fixed vnlRotationMatrix; for (int i = 0; i < 4;++i) for (int j = 0; j < 4;++j) + { transformationMatrix->Element[i][j] = transformation.get(i,j); + if (i < 3 && j < 3) + vnlRotationMatrix.set(i,j,transformation.get(i,j)); + } rotationTransform = itk::ScalableAffineTransform::New(); itk::ScalableAffineTransform::Pointer geometryTransform = img->GetGeometry()->GetIndexToWorldTransform(); img->GetGeometry()->Compose( transformationMatrix); //SetIndexToWorldTransform(geometryTransform); if (dynamic_cast*> (img.GetPointer()) != NULL) { - // apply transformation to image geometry as well as to all gradients !? - } - else - { - // do regular stuff + mitk::DiffusionImage::Pointer diffImages = + dynamic_cast*>(img.GetPointer()); + + mitk::DiffusionImageCorrectionFilter::Pointer correctionFilter = + mitk::DiffusionImageCorrectionFilter::New(); + + correctionFilter->SetImage(diffImages); + // apply transformation to all gradients + correctionFilter->CorrectDirections(vnlRotationMatrix); + //img = correctionFilter->GetOutput(); + } } mitk::BatchedRegistration::TransformType mitk::BatchedRegistration::GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, mitk::Image::Pointer mask) { - + // Handle case that fixed or moving image is a DWI image + mitk::DiffusionImage* fixedDwi = dynamic_cast*> (fixedImage.GetPointer()); + mitk::DiffusionImage* movingDwi = dynamic_cast*> (movingImage.GetPointer()); + itk::B0ImageExtractionImageFilter::Pointer b0Extraction = itk::B0ImageExtractionImageFilter::New(); + if (fixedDwi != NULL) + { + b0Extraction->SetInput(fixedDwi->GetVectorImage()); + b0Extraction->SetDirections(fixedDwi->GetDirections()); + b0Extraction->Update(); + mitk::Image::Pointer tmp = mitk::Image::New(); + tmp->InitializeByItk(b0Extraction->GetOutput()); + tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); + fixedImage = tmp; + } + if (movingDwi != NULL) + { + b0Extraction->SetInput(movingDwi->GetVectorImage()); + b0Extraction->SetDirections(movingDwi->GetDirections()); + b0Extraction->Update(); + mitk::Image::Pointer tmp = mitk::Image::New(); + tmp->InitializeByItk(b0Extraction->GetOutput()); + tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); + movingImage = tmp; + } + // Start registration mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); registrationMethod->SetFixedImage( fixedImage ); if (mask.IsNotNull()) { registrationMethod->SetFixedImageMask(mask); registrationMethod->SetUseFixedImageMask(true); } else { registrationMethod->SetUseFixedImageMask(false); } registrationMethod->SetTransformToRigid(); registrationMethod->SetCrossModalityOn(); + registrationMethod->SetMovingImage(movingImage); registrationMethod->Update(); TransformType transformation; mitk::PyramidImageRegistrationMethod::TransformMatrixType rotationMatrix; rotationMatrix = registrationMethod->GetLastRotationMatrix().transpose(); double param[6]; registrationMethod->GetParameters(param); // first three: euler angles, last three translation for (unsigned int i = 0; i < 3; i++) { for (unsigned int j = 0; j < 3; ++j) { double value = rotationMatrix.get(i,j); transformation.set(i,j, value); } } transformation.set(0,3,-param[3]); transformation.set(1,3,-param[4]); transformation.set(2,3,-param[5]); transformation.set(3,3,1); return transformation; } diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp index 42112e7b0f..9d7a553c3d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp @@ -1,101 +1,117 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP #define MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP #include "mitkDiffusionImageCorrectionFilter.h" #include #include template< typename TPixelType > mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter() { } template< typename TPixelType > typename mitk::DiffusionImageCorrectionFilter::TransformMatrixType mitk::DiffusionImageCorrectionFilter ::GetRotationComponent(const TransformMatrixType &A) { TransformMatrixType B; B = A * A.transpose(); // get the eigenvalues and eigenvectors typedef double MType; vnl_vector< MType > eigvals; vnl_matrix< MType > eigvecs; vnl_symmetric_eigensystem_compute< MType > ( B, eigvecs, eigvals ); vnl_matrix_fixed< MType, 3, 3 > eigvecs_fixed; eigvecs_fixed.set_columns(0, eigvecs ); TransformMatrixType C; C.set_identity(); vnl_vector_fixed< MType, 3 > eigvals_sqrt; for(unsigned int i=0; i<3; i++) { C(i,i) = std::sqrt( eigvals[i] ); } TransformMatrixType S = vnl_inverse( eigvecs_fixed * C * vnl_inverse( eigvecs_fixed )) * A; return S; } template< typename TPixelType > void mitk::DiffusionImageCorrectionFilter ::CorrectDirections( const TransformsVectorType& transformations) { if( m_SourceImage.IsNull() ) { mitkThrow() << " No diffusion image given! "; } GradientDirectionContainerPointerType directions = m_SourceImage->GetDirections(); GradientDirectionContainerPointerType corrected_directions = GradientDirectionContainerType::New(); unsigned int transformed = 0; for(size_t i=0; i< directions->Size(); i++ ) { // skip b-zero images if( directions->ElementAt(i).one_norm() <= 0.0 ) { corrected_directions->push_back( directions->ElementAt(i) ); continue; } GradientDirectionType corrected = GetRotationComponent( transformations.at(transformed)) * directions->ElementAt(i); // store the corrected direction corrected_directions->push_back( corrected ); transformed++; } // replace the old directions with the corrected ones m_SourceImage->SetDirections( corrected_directions ); } +template< typename TPixelType > +void mitk::DiffusionImageCorrectionFilter +::CorrectDirections( const TransformMatrixType& transformation) +{ + if( m_SourceImage.IsNull() ) + { + mitkThrow() << " No diffusion image given! "; + } + TransformsVectorType transfVec; + for (int i=0; i< m_SourceImage->GetDirections()->Size();i++) + { + transfVec.push_back(transformation); + } + this->CorrectDirections(transfVec); +} + #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h index 159cd5321b..7d6c3b3e59 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h @@ -1,90 +1,97 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGECORRECTIONFILTER_H #define MITKDIFFUSIONIMAGECORRECTIONFILTER_H #include "mitkDiffusionImageSource.h" namespace mitk { /** * @class DiffusionImageCorrectionFilter */ template< typename TPixelType > class DiffusionImageCorrectionFilter : public DiffusionImageSource< TPixelType > { public: /** class macros */ mitkClassMacro( DiffusionImageCorrectionFilter, DiffusionImageSource ) itkSimpleNewMacro(Self) typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef vnl_matrix_fixed< double, 3, 3 > TransformMatrixType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef GradientDirectionContainerType::Pointer GradientDirectionContainerPointerType; typedef std::vector< TransformMatrixType > TransformsVectorType; typedef typename Superclass::OutputType DiffusionImageType; typedef typename DiffusionImageType::Pointer DiffusionImageTypePointer; typedef itk::VectorImage ImageType; /** * @brief Set the mitk image ( a 3d+t image ) which is to be reinterpreted as dw image * @param mitkImage */ void SetImage( DiffusionImageTypePointer input ) { m_SourceImage = input; } /** * @brief Correct each gradient direction according to the given transform * * The size of the input is expected to correspond to the count of gradient images in the image. */ void CorrectDirections( const TransformsVectorType& ); + /** + * @brief Correct all gradient directions according to the given transform + * + * This will apply the same rotation to all directions. + */ + void CorrectDirections( const TransformMatrixType& ); + protected: DiffusionImageCorrectionFilter(); virtual ~DiffusionImageCorrectionFilter() {} /** * @brief Get the rotation component following the Finite Strain * * For a given transformation \f$A\f$ its rotation component is defined as \f$ (AA^{T})^{-1/2}\f$. * * The computation first computes \f$ B = AA^T \f$ and then estimates the square root. Square root of * diagonal matrices is defined as * \f$ S = Q * \sqrt{C} * Q^{-1} \f$ with \f$ C \f$ having the eigenvalues on the diagonal. * */ TransformMatrixType GetRotationComponent( const TransformMatrixType& ); DiffusionImageTypePointer m_SourceImage; }; } #include "mitkDiffusionImageCorrectionFilter.cpp" #endif // MITKDIFFUSIONIMAGECORRECTIONFILTER_H