diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
index 6e7bccd182..f1e55fa6ac 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
@@ -1,220 +1,220 @@
/*===================================================================
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 __itkDiffusionMultiShellQballReconstructionImageFilter_h_
#define __itkDiffusionMultiShellQballReconstructionImageFilter_h_
#include <itkImageToImageFilter.h>
namespace itk{
/** \class DiffusionMultiShellQballReconstructionImageFilter
I. Aganj, C. Lenglet, G. Sapiro, E. Yacoub, K. Ugurbil, and N. Harel, “Reconstruction of the orientation distribution function in single and multiple shell q-ball imaging within constant solid angle,” Magnetic Resonance in Medicine, vol. 64, no. 2, pp. 554–566, 2010.
*/
template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections>
class DiffusionMultiShellQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > >
{
public:
typedef DiffusionMultiShellQballReconstructionImageFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
typedef TReferenceImagePixelType ReferencePixelType;
/** GradientImageType
* (e.g. type short)*/
typedef TGradientImagePixelType GradientPixelType;
/** GradientImageType
* 3D VectorImage containing GradientPixelTypes */
typedef VectorImage< GradientPixelType, 3 > GradientImagesType;
/** ODF PixelType */
typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType;
/** ODF ImageType */
typedef Image< OdfPixelType, 3 > OdfImageType;
/** BzeroImageType */
typedef Image< TOdfPixelType, 3 > BZeroImageType;
/** Container to hold gradient directions of the 'n' DW measurements */
typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType;
typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType;
typedef std::map<unsigned int, std::vector<unsigned int> > BValueMap;
typedef std::map<unsigned int, std::vector<unsigned int> >::iterator BValueMapIteraotr;
typedef std::vector<unsigned int> IndiciesVector;
// --------------------------------------------------------------------------------------------//
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter)
/** Get reference image */
virtual typename Superclass::InputImageType * GetInputImage()
{ return ( static_cast< typename Superclass::InputImageType *>(this->ProcessObject::GetInput(0)) ); }
/** Replaces the Input method.
- * Var vols = mitk::DiffusionImage<type>
+ * Var vols = mitk::Image
* -----------------------------------------------------
- * GradientDirectionContainerType-Input gradientDirectionContainer (e.g. vols->GetDirections)
- * GradientImagesType-Input gradientImage (e.g. vols->GetVectorImage)
- * float-Input bvalue (e.g. vols->GetB_Value) */
+ * GradientDirectionContainerType-Input gradientDirectionContainer (e.g. GradientDirectionsContainerProperty)
+ * GradientImagesType-Input gradientImage (e.g. itkVectorImage)
+ * float-Input bvalue (e.g. ReferenceBValueProperty) */
void SetGradientImage( const GradientDirectionContainerType * gradientDirectionContainer, const GradientImagesType *gradientImage , float bvalue);//, std::vector<bool> listOfUserSelctedBValues );
/** Set a BValue Map (key = bvalue, value = indicies splittet for each shell)
* If the input image containes more than three q-shells
* (e.g. b-Values of 0, 1000, 2000, 3000, 4000, ...).
* For the Analytical-Reconstruction it is needed to set a
* BValue Map containing three shells in an arithmetic series
* (e.g. 0, 1000, 2000, 3000).
*/
inline void SetBValueMap(BValueMap map){this->m_BValueMap = map;}
/** Threshold on the reference image data. The output ODF will be a null
* pdf for pixels in the reference image that have a value less than this
* threshold. */
itkSetMacro( Threshold, ReferencePixelType )
itkGetMacro( Threshold, ReferencePixelType )
itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer )
/** Return non-diffusion weighted images */
itkGetMacro( BZeroImage, typename BZeroImageType::Pointer)
/** Factor for Laplacian-Baltrami smoothing of the SH-coefficients*/
itkSetMacro( Lambda, double )
itkGetMacro( Lambda, double )
protected:
DiffusionMultiShellQballReconstructionImageFilter();
~DiffusionMultiShellQballReconstructionImageFilter() { }
void PrintSelf(std::ostream& os, Indent indent) const;
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads );
private:
enum ReconstructionType
{
Mode_Analytical3Shells,
Mode_NumericalNShells,
Mode_Standard1Shell
};
ReconstructionType m_ReconstructionType;
// Interpolation
bool m_Interpolation_Flag;
vnl_matrix< double > * m_Interpolation_SHT1_inv;
vnl_matrix< double > * m_Interpolation_SHT2_inv;
vnl_matrix< double > * m_Interpolation_SHT3_inv;
vnl_matrix< double > * m_TARGET_SH_shell1;
vnl_matrix< double > * m_TARGET_SH_shell2;
vnl_matrix< double > * m_TARGET_SH_shell3;
unsigned int m_MaxDirections;
vnl_matrix< double > * m_CoeffReconstructionMatrix;
vnl_matrix< double > * m_ODFSphericalHarmonicBasisMatrix;
/** container to hold gradient directions */
GradientDirectionContainerType::Pointer m_GradientDirectionContainer;
/** Number of gradient measurements */
unsigned int m_NumberOfGradientDirections;
/** Number of baseline images */
unsigned int m_NumberOfBaselineImages;
/** Threshold on the reference image data */
ReferencePixelType m_Threshold;
typename BZeroImageType::Pointer m_BZeroImage;
typename CoefficientImageType::Pointer m_CoefficientImage;
float m_BValue;
BValueMap m_BValueMap;
double m_Lambda;
bool m_IsHemisphericalArrangementOfGradientDirections;
bool m_IsArithmeticProgession;
void ComputeReconstructionMatrix(IndiciesVector const & refVector);
void ComputeODFSHBasis();
bool CheckDuplicateDiffusionGradients();
bool CheckForDifferingShellDirections();
IndiciesVector GetAllDirections();
void ComputeSphericalHarmonicsBasis(vnl_matrix<double>* QBallReference, vnl_matrix<double>* SHBasisOutput, int Lorder , vnl_matrix<double>* LaplaciaBaltramiOutput =0 , vnl_vector<int>* SHOrderAssociation =0 , vnl_matrix<double> * SHEigenvalues =0);
void Normalize(OdfPixelType & odf );
void S_S0Normalization( vnl_vector<double> & vec, double b0 = 0 );
void DoubleLogarithm(vnl_vector<double> & vec);
double CalculateThreashold(const double value, const double delta);
void Projection1(vnl_vector<double> & vec, double delta = 0.01);
void Projection2( vnl_vector<double> & E1, vnl_vector<double> & E2, vnl_vector<double> & E3, double delta = 0.01);
void Projection3( vnl_vector<double> & A, vnl_vector<double> & alpha, vnl_vector<double> & beta, double delta = 0.01);
void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread);
void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread);
void NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread);
void GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread);
void ComputeSphericalFromCartesian(vnl_matrix<double> * Q, const IndiciesVector & refShell);
//------------------------- VNL-function ------------------------------------
template<typename CurrentValue, typename WntValue>
vnl_vector< WntValue> element_cast (vnl_vector< CurrentValue> const& v1)
{
vnl_vector<WntValue> result(v1.size());
for(unsigned int i = 0 ; i < v1.size(); i++)
result[i] = static_cast< WntValue>(v1[i]);
return result;
}
template<typename type>
double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 )
{
double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm());
return result ;
}
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp"
#endif
#endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
index 32bccf0800..f465a7c29c 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
@@ -1,117 +1,116 @@
/*===================================================================
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 __mitkTeemDiffusionTensor3DReconstructionImageFilter_h__
#define __mitkTeemDiffusionTensor3DReconstructionImageFilter_h__
#include "mitkImage.h"
#include "mitkTensorImage.h"
-#include "mitkDiffusionImage.h"
#include "itkDiffusionTensor3D.h"
namespace mitk
{
enum TeemTensorEstimationMethods{
TeemTensorEstimationMethodsLLS,
TeemTensorEstimationMethodsNLS,
TeemTensorEstimationMethodsWLS,
TeemTensorEstimationMethodsMLE,
};
template< class DiffusionImagePixelType = short,
class TTensorPixelType=float >
class TeemDiffusionTensor3DReconstructionImageFilter : public itk::Object
{
public:
typedef TTensorPixelType TensorPixelType;
typedef itk::Vector<TensorPixelType,7> VectorType;
typedef itk::Image<VectorType,3> VectorImageType;
typedef itk::DiffusionTensor3D<TensorPixelType> TensorType;
typedef itk::Image<TensorType,3 > ItkTensorImageType;
typedef itk::Vector<TensorPixelType,6> ItkTensorVectorType;
typedef itk::Image<ItkTensorVectorType,3> ItkTensorVectorImageType;
typedef DiffusionImagePixelType DiffusionPixelType;
typedef itk::VectorImage< DiffusionPixelType, 3 > DiffusionImageType;
mitkClassMacro( TeemDiffusionTensor3DReconstructionImageFilter,
itk::Object );
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
itkGetMacro(Input,
- typename DiffusionImage<DiffusionImagePixelType>::Pointer);
+ mitk::Image::Pointer);
itkSetMacro(Input,
- typename DiffusionImage<DiffusionImagePixelType>::Pointer);
+ mitk::Image::Pointer);
itkGetMacro(EstimateErrorImage, bool);
itkSetMacro(EstimateErrorImage, bool);
itkGetMacro(Sigma, float);
itkSetMacro(Sigma, float);
itkGetMacro(EstimationMethod, TeemTensorEstimationMethods);
itkSetMacro(EstimationMethod, TeemTensorEstimationMethods);
itkGetMacro(NumIterations, int);
itkSetMacro(NumIterations, int);
itkGetMacro(ConfidenceThreshold, double);
itkSetMacro(ConfidenceThreshold, double);
itkGetMacro(ConfidenceFuzzyness, float);
itkSetMacro(ConfidenceFuzzyness, float);
itkGetMacro(MinPlausibleValue, double);
itkSetMacro(MinPlausibleValue, double);
itkGetMacro(Output, mitk::TensorImage::Pointer);
itkGetMacro(OutputItk, mitk::TensorImage::Pointer);
// do the work
virtual void Update();
protected:
TeemDiffusionTensor3DReconstructionImageFilter();
virtual ~TeemDiffusionTensor3DReconstructionImageFilter();
- typename DiffusionImage<DiffusionImagePixelType>::Pointer m_Input;
+ mitk::Image::Pointer m_Input;
bool m_EstimateErrorImage;
float m_Sigma;
TeemTensorEstimationMethods m_EstimationMethod;
int m_NumIterations;
double m_ConfidenceThreshold;
float m_ConfidenceFuzzyness;
double m_MinPlausibleValue;
mitk::TensorImage::Pointer m_Output;
mitk::TensorImage::Pointer m_OutputItk;
mitk::Image::Pointer m_ErrorImage;
};
}
#include "mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp"
#endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp
index d11c505925..3f6fbc61a9 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp
@@ -1,327 +1,316 @@
/*===================================================================
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 MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP
-#define MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP
+#ifndef MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
+#define MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
#include "mitkDWIHeadMotionCorrectionFilter.h"
#include "itkSplitDWImageFilter.h"
#include "itkB0ImageExtractionToSeparateImageFilter.h"
#include "mitkImageTimeSelector.h"
#include "mitkPyramidImageRegistrationMethod.h"
//#include "mitkRegistrationMethodITK4.h"
-#include "mitkImageToDiffusionImageSource.h"
+#include <mitkDiffusionPropertyHelper.h>
#include "mitkDiffusionImageCorrectionFilter.h"
#include <mitkImageWriteAccessor.h>
+#include <mitkProperties.h>
#include <vector>
#include "mitkIOUtil.h"
#include <itkImage.h>
-template< typename DiffusionPixelType>
-mitk::DWIHeadMotionCorrectionFilter<DiffusionPixelType>
- ::DWIHeadMotionCorrectionFilter()
+mitk::DWIHeadMotionCorrectionFilter::DWIHeadMotionCorrectionFilter()
: m_CurrentStep(0),
m_Steps(100),
m_IsInValidState(true),
m_AbortRegistration(false),
m_AverageUnweighted(true)
{
}
-template< typename DiffusionPixelType>
-void mitk::DWIHeadMotionCorrectionFilter<DiffusionPixelType>
-::GenerateData()
+
+void mitk::DWIHeadMotionCorrectionFilter::GenerateData()
{
typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType;
- DiffusionImageType* input = const_cast<DiffusionImageType*>(this->GetInput(0));
+ InputImageType* input = const_cast<InputImageType*>(this->GetInput(0));
+
+ ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New();
+ mitk::CastToItkImage(input, itkVectorImagePointer);
typedef mitk::PyramidImageRegistrationMethod RegistrationMethod;
// typedef mitk::RegistrationMethodITKV4 RegistrationMethod
-
- unsigned int numberOfSteps = input->GetVectorImage()->GetNumberOfComponentsPerPixel () ;
+ unsigned int numberOfSteps = itkVectorImagePointer->GetNumberOfComponentsPerPixel () ;
m_Steps = numberOfSteps;
//
// (1) Extract the b-zero images to a 3d+t image, register them by NCorr metric and
// rigid registration : they will then be used are reference image for registering
// the gradient images
//
typedef itk::B0ImageExtractionToSeparateImageFilter< DiffusionPixelType, DiffusionPixelType> B0ExtractorType;
- typename B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New();
- b0_extractor->SetInput( input->GetVectorImage() );
- b0_extractor->SetDirections( input->GetDirections() );
+ B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New();
+ b0_extractor->SetInput( itkVectorImagePointer );
+ b0_extractor->SetDirections( static_cast<mitk::GradientDirectionsProperty*>( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() );
b0_extractor->Update();
mitk::Image::Pointer b0Image = mitk::Image::New();
b0Image->InitializeByItk( b0_extractor->GetOutput() );
b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(),
mitk::Image::CopyMemory );
// (2.1) Use the extractor to access the extracted b0 volumes
mitk::ImageTimeSelector::Pointer t_selector =
mitk::ImageTimeSelector::New();
t_selector->SetInput( b0Image );
t_selector->SetTimeNr(0);
t_selector->Update();
// first unweighted image as reference space for the registration
mitk::Image::Pointer b0referenceImage = t_selector->GetOutput();
const unsigned int numberOfb0Images = b0Image->GetTimeSteps();
// register b-zeros only if the flag to average is set, use the first one otherwise
if( m_AverageUnweighted && numberOfb0Images > 1)
{
RegistrationMethod::Pointer registrationMethod = RegistrationMethod::New();
registrationMethod->SetFixedImage( b0referenceImage );
registrationMethod->SetTransformToRigid();
// the unweighted images are of same modality
registrationMethod->SetCrossModalityOff();
// use the advanced (windowed sinc) interpolation
registrationMethod->SetUseAdvancedInterpolation(true);
// Initialize the temporary output image
mitk::Image::Pointer registeredB0Image = b0Image->Clone();
mitk::ImageTimeSelector::Pointer t_selector2 =
mitk::ImageTimeSelector::New();
t_selector2->SetInput( b0Image );
for( unsigned int i=1; i<numberOfb0Images; i++)
{
m_CurrentStep = i + 1;
if( m_AbortRegistration == true) {
m_IsInValidState = false;
mitkThrow() << "Stopped by user.";
};
t_selector2->SetTimeNr(i);
t_selector2->Update();
registrationMethod->SetMovingImage( t_selector2->GetOutput() );
try
{
MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration";
registrationMethod->Update();
}
catch( const itk::ExceptionObject& e)
{
m_IsInValidState = false;
mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what();
}
// import volume to the inter-results
mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage());
registeredB0Image->SetImportVolume( imac.GetData(),
i, 0, mitk::Image::ReferenceMemory );
}
// use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file
AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage );
}
//
// (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images
//
- typename SplitFilterType::Pointer split_filter = SplitFilterType::New();
- split_filter->SetInput (input->GetVectorImage() );
- split_filter->SetExtractAllAboveThreshold(8, input->GetBValueMap() );
+ SplitFilterType::Pointer split_filter = SplitFilterType::New();
+ split_filter->SetInput (itkVectorImagePointer );
+ split_filter->SetExtractAllAboveThreshold(8, static_cast<mitk::BValueMapProperty*>(input->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap() );
try
{
split_filter->Update();
}
catch( const itk::ExceptionObject &e)
{
m_IsInValidState = false;
mitkThrow() << " Caught exception from SplitImageFilter : " << e.what();
}
mitk::Image::Pointer splittedImage = mitk::Image::New();
splittedImage->InitializeByItk( split_filter->GetOutput() );
splittedImage->SetImportChannel( split_filter->GetOutput()->GetBufferPointer(),
mitk::Image::CopyMemory );
//
// (3) Use again the time-selector to access the components separately in order
// to perform the registration of Image -> unweighted reference
//
RegistrationMethod::Pointer weightedRegistrationMethod
= RegistrationMethod::New();
weightedRegistrationMethod->SetTransformToAffine();
weightedRegistrationMethod->SetCrossModalityOn();
//
// - (3.1) Set the reference image
// - a single b0 image
// - average over the registered b0 images if multiple present
//
weightedRegistrationMethod->SetFixedImage( b0referenceImage );
// use the advanced (windowed sinc) interpolation
weightedRegistrationMethod->SetUseAdvancedInterpolation(true);
weightedRegistrationMethod->SetVerboseOn();
//
// - (3.2) Register all timesteps in the splitted image onto the first reference
//
unsigned int maxImageIdx = splittedImage->GetTimeSteps();
mitk::TimeGeometry* tsg = splittedImage->GetTimeGeometry();
mitk::ProportionalTimeGeometry* ptg = dynamic_cast<ProportionalTimeGeometry*>(tsg);
ptg->Expand(maxImageIdx+1);
ptg->SetTimeStepGeometry( ptg->GetGeometryForTimeStep(0), maxImageIdx );
mitk::Image::Pointer registeredWeighted = mitk::Image::New();
registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg );
// insert the first unweighted reference as the first volume
// in own scope to release the accessor asap after copy
{
mitk::ImageWriteAccessor imac(b0referenceImage);
registeredWeighted->SetImportVolume( imac.GetData(),
0,0, mitk::Image::CopyMemory );
}
// mitk::Image::Pointer registeredWeighted = splittedImage->Clone();
// this time start at 0, we have only gradient images in the 3d+t file
// the reference image comes form an other image
mitk::ImageTimeSelector::Pointer t_selector_w =
mitk::ImageTimeSelector::New();
t_selector_w->SetInput( splittedImage );
// store the rotation parts of the transformations in a vector
typedef RegistrationMethod::TransformMatrixType MatrixType;
std::vector< MatrixType > estimated_transforms;
for( unsigned int i=0; i<maxImageIdx; i++)
{
m_CurrentStep = numberOfb0Images + i + 1;
if( m_AbortRegistration == true) {
m_IsInValidState = false;
mitkThrow() << "Stopped by user.";
};
t_selector_w->SetTimeNr(i);
t_selector_w->Update();
weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() );
try
{
MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration";
weightedRegistrationMethod->Update();
}
catch( const itk::ExceptionObject& e)
{
m_IsInValidState = false;
mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what();
}
// allow expansion
mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage());
registeredWeighted->SetImportVolume( imac.GetData(),
i+1, 0, mitk::Image::CopyMemory);
estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() );
}
//
// (4) Cast the resulting image back to an diffusion weighted image
//
- typename DiffusionImageType::GradientDirectionContainerType *gradients = input->GetDirections();
- typename DiffusionImageType::GradientDirectionContainerType::Pointer gradients_new =
- DiffusionImageType::GradientDirectionContainerType::New();
- typename DiffusionImageType::GradientDirectionType bzero_vector;
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType *gradients = static_cast<mitk::GradientDirectionsProperty*>( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradients_new =
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New();
+ mitk::DiffusionPropertyHelper::GradientDirectionType bzero_vector;
bzero_vector.fill(0);
// compose the direction vector
// - no direction for the first image
// - correct ordering of the directions based on the index list
gradients_new->push_back( bzero_vector );
- typename SplitFilterType::IndexListType index_list = split_filter->GetIndexList();
- typename SplitFilterType::IndexListType::const_iterator lIter = index_list.begin();
+ SplitFilterType::IndexListType index_list = split_filter->GetIndexList();
+ SplitFilterType::IndexListType::const_iterator lIter = index_list.begin();
while( lIter != index_list.end() )
{
gradients_new->push_back( gradients->at( *lIter ) );
++lIter;
}
- typename mitk::ImageToDiffusionImageSource< DiffusionPixelType >::Pointer caster =
- mitk::ImageToDiffusionImageSource< DiffusionPixelType >::New();
-
- caster->SetImage( registeredWeighted );
- caster->SetBValue( input->GetReferenceBValue() );
- caster->SetGradientDirections( gradients_new.GetPointer() );
-
- try
- {
- caster->Update();
- }
- catch( const itk::ExceptionObject& e)
- {
- m_IsInValidState = false;
- MITK_ERROR << "Casting back to diffusion image failed: ";
- mitkThrow() << "Subprocess failed with exception: " << e.what();
- }
+ registeredWeighted->SetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradients_new.GetPointer() ));
+ registeredWeighted->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast<mitk::FloatProperty*>(input->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) );
//
// (5) Adapt the gradient directions according to the estimated transforms
//
- typedef mitk::DiffusionImageCorrectionFilter< DiffusionPixelType > CorrectionFilterType;
- typename CorrectionFilterType::Pointer corrector = CorrectionFilterType::New();
+ typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType;
+ CorrectionFilterType::Pointer corrector = CorrectionFilterType::New();
- OutputImagePointerType output = caster->GetOutput();
+ mitk::Image::Pointer output = registeredWeighted;
corrector->SetImage( output );
corrector->CorrectDirections( estimated_transforms );
//
// (6) Pass the corrected image to the filters output port
//
m_CurrentStep += 1;
- this->GetOutput()->SetVectorImage(output->GetVectorImage());
- this->GetOutput()->SetReferenceBValue(output->GetReferenceBValue());
- this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame());
- this->GetOutput()->SetDirections(output->GetDirections());
- this->GetOutput()->InitializeFromVectorImage();
+
+ this->GetOutput()->Initialize( output );
+
+ this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast<mitk::GradientDirectionsProperty*>( output->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) );
+ this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast<mitk::MeasurementFrameProperty*>( output->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) );
+ this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast<mitk::FloatProperty*>(output->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) );
+
+ mitk::DiffusionPropertyHelper propertyHelper( this->GetOutput() );
+ propertyHelper.InitializeImage();
this->GetOutput()->Modified();
}
-#endif // MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP
+#endif // MITKDWIHEADMOTIONCORRECTIONFILTER_CPP
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h
index d8065ff27d..e159597b08 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h
@@ -1,132 +1,131 @@
/*===================================================================
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 MITKDWIHEADMOTIONCORRECTIONFILTER_H
#define MITKDWIHEADMOTIONCORRECTIONFILTER_H
-#include "mitkDiffusionImageToDiffusionImageFilter.h"
+#include "mitkImageToImageFilter.h"
#include <itkAccumulateImageFilter.h>
#include <itkExtractImageFilter.h>
#include "mitkITKImageImport.h"
+#include <itkVectorImage.h>
+#include <MitkDiffusionCoreExports.h>
namespace mitk
{
/**
* @class DWIHeadMotionCorrectionFilter
*
* @brief Performs standard head-motion correction by using affine registration of the gradient images.
*
* (Head) motion correction is a essential pre-processing step before performing any further analysis of a diffusion-weighted
* images since all model fits ( tensor, QBI ) rely on an aligned diffusion-weighted dataset. The correction is done in two steps. First the
* unweighted images ( if multiple present ) are separately registered on the first one by means of rigid registration and normalized correlation
* as error metric. Second, the weighted gradient images are registered to the unweighted reference ( computed as average from the aligned images from first step )
* by an affine transformation using the MattesMutualInformation metric as optimizer guidance.
*
*/
-template< typename DiffusionPixelType>
-class DWIHeadMotionCorrectionFilter
- : public DiffusionImageToDiffusionImageFilter< DiffusionPixelType >
+
+class MitkDiffusionCore_EXPORT DWIHeadMotionCorrectionFilter
+ : public ImageToImageFilter
{
public:
// class macros
mitkClassMacro( DWIHeadMotionCorrectionFilter,
- DiffusionImageToDiffusionImageFilter<DiffusionPixelType> )
+ ImageToImageFilter )
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
itkGetMacro( Steps, unsigned long )
itkGetMacro( CurrentStep, unsigned long )
itkGetMacro( IsInValidState, bool)
itkSetMacro( AbortRegistration, bool )
itkSetMacro( AverageUnweighted, bool )
// public typedefs
- typedef typename Superclass::InputImageType DiffusionImageType;
- typedef typename Superclass::InputImagePointerType DiffusionImagePointerType;
-
- typedef typename Superclass::OutputImageType OutputImageType;
- typedef typename Superclass::OutputImagePointerType OutputImagePointerType;
+ typedef short DiffusionPixelType;
+ typedef Superclass::InputImageType InputImageType;
+ typedef Superclass::OutputImageType OutputImageType;
+ typedef itk::VectorImage<DiffusionPixelType, 3> ITKDiffusionImageType;
protected:
DWIHeadMotionCorrectionFilter();
virtual ~DWIHeadMotionCorrectionFilter() {}
virtual void GenerateData();
unsigned long m_CurrentStep;
unsigned long m_Steps;
bool m_IsInValidState; ///< Whether the filter is in a valid state, false if error occured
bool m_AbortRegistration; ///< set flag to abort
bool m_AverageUnweighted;
};
/**
* @brief Averages an 3d+t image along the time axis.
*
* The method uses the AccumulateImageFilter as provided by ITK and collapses the given 3d+t image
* to an 3d image while computing the average over the time axis for each of the spatial voxels.
*/
template< typename TPixel, unsigned int VDimensions>
static void ItkAccumulateFilter(
const itk::Image< TPixel, VDimensions>* image,
mitk::Image::Pointer& output)
{
// input 3d+t --> output 3d
typedef itk::Image< TPixel, 4> InputItkType;
typedef itk::Image< TPixel, 3> OutputItkType;
// the accumulate filter requires the same dimension in output and input image
typedef typename itk::AccumulateImageFilter< InputItkType, InputItkType > FilterType;
typename FilterType::Pointer filter = FilterType::New();
filter->SetInput( image );
filter->SetAccumulateDimension( 3 );
filter->SetAverage( true );
// we need to extract the volume to reduce the size from 4 to 3 for further processing
typedef typename itk::ExtractImageFilter< InputItkType, OutputItkType > ExtractFilterType;
typename ExtractFilterType::Pointer extractor = ExtractFilterType::New();
extractor->SetInput( filter->GetOutput() );
extractor->SetDirectionCollapseToIdentity();
typename InputItkType::RegionType extractRegion = image->GetLargestPossibleRegion();
// crop out the time axis
extractRegion.SetSize( 3, 0);
extractor->SetExtractionRegion( extractRegion );
try
{
extractor->Update();
}
catch( const itk::ExceptionObject& e)
{
mitkThrow() << " Exception while averaging: " << e.what();
}
output = mitk::GrabItkImageMemory( extractor->GetOutput() );
}
} //end namespace mitk
-#include "mitkDWIHeadMotionCorrectionFilter.cpp"
-
#endif // MITKDWIHEADMOTIONCORRECTIONFILTER_H
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp
index 0e1982705b..2f488f4334 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp
@@ -1,196 +1,216 @@
#include "mitkRegistrationWrapper.h"
#include "mitkPyramidImageRegistrationMethod.h"
-#include "mitkDiffusionImage.h"
+#include "mitkImage.h"
#include <mitkDiffusionImageCorrectionFilter.h>
#include "itkB0ImageExtractionImageFilter.h"
#include <itkResampleImageFilter.h>
#include <itkWindowedSincInterpolateImageFunction.h>
#include <itkLinearInterpolateImageFunction.h>
#include <itkNearestNeighborExtrapolateImageFunction.h>
#include <itkBSplineInterpolateImageFunction.h>
+#include <mitkDiffusionPropertyHelper.h>
#include <vnl/vnl_inverse.h>
void mitk::RegistrationWrapper::ApplyTransformationToImage(mitk::Image::Pointer img, const mitk::RegistrationWrapper::RidgidTransformType &transformation,double* offset, mitk::Image* resampleReference, bool binary)
{
- typedef mitk::DiffusionImage<short> DiffusionImageType;
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientDirections =
+ static_cast<mitk::GradientDirectionsProperty*>( img->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
- if (dynamic_cast<DiffusionImageType*> (img.GetPointer()) == NULL)
+ if (gradientDirections.IsNull() || ( gradientDirections->Size() == 0) )
{
ItkImageType::Pointer itkImage = ItkImageType::New();
CastToItkImage(img, itkImage);
typedef itk::Euler3DTransform< double > RigidTransformType;
RigidTransformType::Pointer rtransform = RigidTransformType::New();
RigidTransformType::ParametersType parameters(RigidTransformType::ParametersDimension);
for (int i = 0; i<6;++i)
parameters[i] = transformation[i];
rtransform->SetParameters( parameters );
mitk::Point3D origin = itkImage->GetOrigin();
origin[0]-=offset[0];
origin[1]-=offset[1];
origin[2]-=offset[2];
mitk::Point3D newOrigin = rtransform->GetInverseTransform()->TransformPoint(origin);
itk::Matrix<double,3,3> dir = itkImage->GetDirection();
itk::Matrix<double,3,3> transM ( vnl_inverse(rtransform->GetMatrix().GetVnlMatrix()));
itk::Matrix<double,3,3> newDirection = transM * dir;
itkImage->SetOrigin(newOrigin);
itkImage->SetDirection(newDirection);
// Perform Resampling if reference image is provided
if (resampleReference != NULL)
{
typedef itk::ResampleImageFilter<ItkImageType, ItkImageType> ResampleFilterType;
ItkImageType::Pointer itkReference = ItkImageType::New();
CastToItkImage(resampleReference,itkReference);
typedef itk::Function::WelchWindowFunction<4> WelchWindowFunction;
typedef itk::WindowedSincInterpolateImageFunction< ItkImageType, 4,WelchWindowFunction> WindowedSincInterpolatorType;
WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New();
typedef itk::LinearInterpolateImageFunction< ItkImageType> LinearInterpolatorType;
LinearInterpolatorType::Pointer lin_interpolator = LinearInterpolatorType::New();
typedef itk::NearestNeighborInterpolateImageFunction< ItkImageType, double > NearestNeighborInterpolatorType;
NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New();
typedef itk::BSplineInterpolateImageFunction< ItkImageType, double > BSplineInterpolatorType;
BSplineInterpolatorType::Pointer bSpline_interpolator = BSplineInterpolatorType::New();
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetInput(itkImage);
resampler->SetReferenceImage( itkReference );
resampler->UseReferenceImageOn();
if (binary)
resampler->SetInterpolator(nn_interpolator);
else
resampler->SetInterpolator(lin_interpolator);
resampler->Update();
GrabItkImageMemory(resampler->GetOutput(), img);
}
else
{
// !! CastToItk behaves very differently depending on the original data type
// if the target type is the same as the original, only a pointer to the data is set
// and an additional GrabItkImageMemory will cause a segfault when the image is destroyed
// GrabItkImageMemory - is not necessary in this case since we worked on the original data
// See Bug 17538.
if (img->GetPixelType().GetComponentTypeAsString() != "double")
img = GrabItkImageMemory(itkImage);
}
}
else
{
- DiffusionImageType::Pointer diffImages = dynamic_cast<DiffusionImageType*>(img.GetPointer());
+ mitk::Image::Pointer diffImages = dynamic_cast<mitk::Image*>(img.GetPointer());
typedef itk::Euler3DTransform< double > RigidTransformType;
RigidTransformType::Pointer rtransform = RigidTransformType::New();
RigidTransformType::ParametersType parameters(RigidTransformType::ParametersDimension);
for (int i = 0; i<6;++i)
{
parameters[i] = transformation[i];
}
rtransform->SetParameters( parameters );
- mitk::Point3D b0origin = diffImages->GetVectorImage()->GetOrigin();
+ typedef itk::VectorImage<short, 3> ITKDiffusionImageType;
+ ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New();
+ mitk::CastToItkImage(diffImages, itkVectorImagePointer);
+
+ mitk::Point3D b0origin = itkVectorImagePointer->GetOrigin();
b0origin[0]-=offset[0];
b0origin[1]-=offset[1];
b0origin[2]-=offset[2];
mitk::Point3D newOrigin = rtransform->GetInverseTransform()->TransformPoint(b0origin);
- itk::Matrix<double,3,3> dir = diffImages->GetVectorImage()->GetDirection();
+ itk::Matrix<double,3,3> dir = itkVectorImagePointer->GetDirection();
itk::Matrix<double,3,3> transM ( vnl_inverse(rtransform->GetMatrix().GetVnlMatrix()));
itk::Matrix<double,3,3> newDirection = transM * dir;
- diffImages->GetVectorImage()->SetOrigin(newOrigin);
- diffImages->GetVectorImage()->SetDirection(newDirection);
+ itkVectorImagePointer->SetOrigin(newOrigin);
+ itkVectorImagePointer->SetDirection(newDirection);
diffImages->Modified();
- mitk::DiffusionImageCorrectionFilter<short>::Pointer correctionFilter = mitk::DiffusionImageCorrectionFilter<short>::New();
+ mitk::DiffusionImageCorrectionFilter::Pointer correctionFilter = mitk::DiffusionImageCorrectionFilter::New();
// For Diff. Images: Need to rotate the gradients (works in-place)
correctionFilter->SetImage(diffImages);
correctionFilter->CorrectDirections(transM.GetVnlMatrix());
img = diffImages;
}
}
void mitk::RegistrationWrapper::GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, RidgidTransformType transformation,double* offset, bool useSameOrigin, mitk::Image* mask)
{
// Handle the case that fixed/moving image is a DWI image
- mitk::DiffusionImage<short>* fixedDwi = dynamic_cast<mitk::DiffusionImage<short>*> (fixedImage.GetPointer());
- mitk::DiffusionImage<short>* movingDwi = dynamic_cast<mitk::DiffusionImage<short>*> (movingImage.GetPointer());
+ mitk::Image* fixedDwi = dynamic_cast<mitk::Image*> (fixedImage.GetPointer());
+ mitk::Image* movingDwi = dynamic_cast<mitk::Image*> (movingImage.GetPointer());
itk::B0ImageExtractionImageFilter<short,short >::Pointer b0Extraction = itk::B0ImageExtractionImageFilter<short,short>::New();
offset[0]=offset[1]=offset[2]=0;
- if (fixedDwi != NULL)
+
+ typedef itk::VectorImage<short, 3> ITKDiffusionImageType;
+ ITKDiffusionImageType::Pointer itkFixedDwiPointer = ITKDiffusionImageType::New();
+ mitk::CastToItkImage(fixedDwi, itkFixedDwiPointer);
+
+ ITKDiffusionImageType::Pointer itkMovingDwiPointer = ITKDiffusionImageType::New();
+ mitk::CastToItkImage(movingDwi, itkMovingDwiPointer);
+
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer fixedGradientDirections =
+ static_cast<mitk::GradientDirectionsProperty*>( fixedDwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
+
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer movingGradientDirections =
+ static_cast<mitk::GradientDirectionsProperty*>( movingDwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
+
+ if (fixedGradientDirections.IsNull() || ( fixedGradientDirections->Size() == 0))
{ // Set b0 extraction as fixed image
- b0Extraction->SetInput(fixedDwi->GetVectorImage());
- b0Extraction->SetDirections(fixedDwi->GetDirections());
+ b0Extraction->SetInput( itkFixedDwiPointer );
+ b0Extraction->SetDirections(fixedGradientDirections);
b0Extraction->Update();
mitk::Image::Pointer tmp = mitk::Image::New();
tmp->InitializeByItk(b0Extraction->GetOutput());
tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer());
fixedImage = tmp;
}
- if (movingDwi != NULL)
+ if (movingGradientDirections || ( movingGradientDirections->Size() == 0))
{ // Set b0 extraction as moving image
- b0Extraction->SetInput(movingDwi->GetVectorImage());
- b0Extraction->SetDirections(movingDwi->GetDirections());
+ b0Extraction->SetInput( itkMovingDwiPointer );
+ b0Extraction->SetDirections(movingGradientDirections);
b0Extraction->Update();
mitk::Image::Pointer tmp = mitk::Image::New();
tmp->InitializeByItk(b0Extraction->GetOutput());
tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer());
movingImage = tmp;
}
// align the offsets of the two images. this is done to avoid non-overlapping initialization
Point3D origin = fixedImage->GetGeometry()->GetOrigin();
Point3D originMoving = movingImage->GetGeometry()->GetOrigin();
mitk::Image::Pointer tmpImage = movingImage->Clone();
if (useSameOrigin)
{
offset[0] = originMoving[0]-origin[0];
offset[1] = originMoving[1]-origin[1];
offset[2] = originMoving[2]-origin[2];
tmpImage->GetGeometry()->SetOrigin(origin);
}
// Start registration
mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New();
registrationMethod->SetFixedImage( fixedImage );
if (mask != NULL)
{
registrationMethod->SetFixedImageMask(mask);
registrationMethod->SetUseFixedImageMask(true);
}
else
{
registrationMethod->SetUseFixedImageMask(false);
}
registrationMethod->SetTransformToRigid();
registrationMethod->SetCrossModalityOn();
registrationMethod->SetMovingImage(tmpImage);
registrationMethod->Update();
registrationMethod->GetParameters(transformation); // first three: euler angles, last three translation
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h
index deead178e9..7b7802c3cc 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h
@@ -1,80 +1,80 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
/*===================================================================
This file is based heavily on a corresponding ITK filter.
===================================================================*/
#ifndef __itkAdcImageFilter_h_
#define __itkAdcImageFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
-#include <mitkDiffusionImage.h>
+#include <mitkDiffusionPropertyHelper.h>
namespace itk{
/** \class AdcImageFilter
*/
template< class TInPixelType, class TOutPixelType >
class AdcImageFilter :
public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > >
{
public:
typedef AdcImageFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > Superclass;
- typedef mitk::DiffusionImage< short >::GradientDirectionType GradientDirectionType;
- typedef mitk::DiffusionImage< short >::GradientDirectionContainerType::Pointer GradientContainerType;
+ typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType;
+ typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(AdcImageFilter, ImageToImageFilter)
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
itkSetMacro( B_value, double )
itkSetMacro( GradientDirections, GradientContainerType )
protected:
AdcImageFilter();
~AdcImageFilter() {}
void PrintSelf(std::ostream& os, Indent indent) const;
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType);
double m_B_value;
GradientContainerType m_GradientDirections;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkAdcImageFilter.txx"
#endif
#endif //__itkAdcImageFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
index cdb43eec2e..c268bb5e40 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h
@@ -1,100 +1,100 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
/*===================================================================
This file is based heavily on a corresponding ITK filter.
===================================================================*/
#ifndef __itkDwiNormilzationFilter_h_
#define __itkDwiNormilzationFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
-#include <mitkDiffusionImage.h>
#include <itkLabelStatisticsImageFilter.h>
#include <itkImageRegionIteratorWithIndex.h>
#include <itkShiftScaleImageFilter.h>
+#include <mitkDiffusionPropertyHelper.h>
namespace itk{
/** \class DwiNormilzationFilter
* \brief Normalizes the data vectors either using the specified reference value or the voxelwise baseline value.
*/
template< class TInPixelType >
class DwiNormilzationFilter :
public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > >
{
public:
typedef DwiNormilzationFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(DwiNormilzationFilter, ImageToImageFilter)
typedef itk::Image< double, 3 > DoubleImageType;
typedef itk::Image< unsigned char, 3 > UcharImageType;
typedef itk::Image< unsigned short, 3 > BinImageType;
typedef itk::Image< TInPixelType, 3 > TInPixelImageType;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
- typedef mitk::DiffusionImage< short >::GradientDirectionType GradientDirectionType;
- typedef mitk::DiffusionImage< short >::GradientDirectionContainerType::Pointer GradientContainerType;
+ typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType;
+ typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType;
typedef itk::LabelStatisticsImageFilter< TInPixelImageType,BinImageType > StatisticsFilterType;
typedef itk::Statistics::Histogram< typename TInPixelImageType::PixelType > HistogramType;
typedef typename HistogramType::MeasurementType MeasurementType;
typedef itk::ShiftScaleImageFilter<TInPixelImageType, DoubleImageType> ShiftScaleImageFilterType;
itkSetMacro( GradientDirections, GradientContainerType )
itkSetMacro( ScalingFactor, TInPixelType )
itkSetMacro( UseGlobalReference, bool )
itkSetMacro( MaskImage, UcharImageType::Pointer )
itkSetMacro( Reference, double )
protected:
DwiNormilzationFilter();
~DwiNormilzationFilter() {}
void PrintSelf(std::ostream& os, Indent indent) const;
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType);
UcharImageType::Pointer m_MaskImage;
GradientContainerType m_GradientDirections;
int m_B0Index;
TInPixelType m_ScalingFactor;
bool m_UseGlobalReference;
double m_Reference;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkDwiNormilzationFilter.txx"
#endif
#endif //__itkDwiNormilzationFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h
index dd75de3f14..9c80bb27e4 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h
@@ -1,123 +1,126 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_
#define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_
#include <itkImageToImageFilter.h>
#include <itkVectorImage.h>
#include <itkPointShell.h>
namespace itk
{
/**
* \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */
template <class TInputScalarType, class TOutputScalarType>
class ElectrostaticRepulsionDiffusionGradientReductionFilter
: public ImageToImageFilter<itk::VectorImage<TInputScalarType,3>, itk::VectorImage<TOutputScalarType,3> >
{
public:
typedef ElectrostaticRepulsionDiffusionGradientReductionFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< itk::VectorImage<TOutputScalarType,3>, itk::VectorImage<TOutputScalarType,3> >
Superclass;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(ElectrostaticRepulsionDiffusionGradientReductionFilter, ImageToImageFilter)
typedef TInputScalarType InputScalarType;
typedef itk::VectorImage<InputScalarType,3> InputImageType;
typedef typename InputImageType::PixelType InputPixelType;
typedef TOutputScalarType OutputScalarType;
typedef itk::VectorImage<OutputScalarType,3> OutputImageType;
typedef typename OutputImageType::PixelType OutputPixelType;
typedef OutputScalarType BaselineScalarType;
typedef BaselineScalarType BaselinePixelType;
typedef typename itk::Image<BaselinePixelType,3> BaselineImageType;
typedef vnl_vector_fixed< double, 3 > GradientDirectionType;
typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType;
typedef std::vector<unsigned int> IndicesVector;
typedef std::map<unsigned int, IndicesVector> BValueMap;
itkGetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer)
itkSetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer)
itkGetMacro(GradientDirections, GradientDirectionContainerType::Pointer)
itkSetMacro(GradientDirections, GradientDirectionContainerType::Pointer)
IndicesVector GetUsedGradientIndices(){return m_UsedGradientIndices;}
void SetOriginalBValueMap(BValueMap inp){m_OriginalBValueMap = inp;}
void SetShellSelectionBValueMap(BValueMap inp){m_InputBValueMap = inp;}
void SetNumGradientDirections(std::vector<unsigned int> numDirs){m_NumGradientDirections = numDirs;}
+
+ void UpdateOutputInformation();
+
protected:
ElectrostaticRepulsionDiffusionGradientReductionFilter();
~ElectrostaticRepulsionDiffusionGradientReductionFilter() {}
void GenerateData();
double Costs(); ///< calculates electrostatic energy of current direction set
GradientDirectionContainerType::Pointer m_GradientDirections; ///< container for the subsampled output gradient directions
GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions
IndicesVector m_UsedGradientIndices;
IndicesVector m_UnusedGradientIndices;
IndicesVector m_BaselineImageIndices;
BValueMap m_OriginalBValueMap;
BValueMap m_InputBValueMap;
std::vector<unsigned int> m_NumGradientDirections;
};
} // end of namespace
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx"
#endif
#endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx
index 3b601cdb0a..7253e35b42 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx
@@ -1,247 +1,261 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_
#define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_
#endif
#define _USE_MATH_DEFINES
#include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.h"
#include <math.h>
#include <time.h>
#include <itkImageRegionIterator.h>
#include <itkImageRegion.h>
namespace itk
{
template <class TInputScalarType, class TOutputScalarType>
ElectrostaticRepulsionDiffusionGradientReductionFilter<TInputScalarType, TOutputScalarType>
::ElectrostaticRepulsionDiffusionGradientReductionFilter()
{
this->SetNumberOfRequiredInputs( 1 );
}
template <class TInputScalarType, class TOutputScalarType>
double
ElectrostaticRepulsionDiffusionGradientReductionFilter<TInputScalarType, TOutputScalarType>
::Costs()
{
double costs = 2*M_PI;
for (IndicesVector::iterator it = m_UsedGradientIndices.begin(); it!=m_UsedGradientIndices.end(); ++it)
{
for (IndicesVector::iterator it2 = m_UsedGradientIndices.begin(); it2!=m_UsedGradientIndices.end(); ++it2)
if (it != it2)
{
vnl_vector_fixed<double,3> v1 = m_OriginalGradientDirections->at(*it);
vnl_vector_fixed<double,3> v2 = m_OriginalGradientDirections->at(*it2);
v1.normalize(); v2.normalize();
double temp = dot_product(v1,v2);
if (temp>1)
temp = 1;
else if (temp<-1)
temp = -1;
double angle = acos(temp);
if (angle<costs)
costs = angle;
temp = dot_product(-v1,v2);
if (temp>1)
temp = 1;
else if (temp<-1)
temp = -1;
angle = acos(temp);
if (angle<costs)
costs = angle;
}
}
return costs;
}
template <class TInputScalarType, class TOutputScalarType>
void
ElectrostaticRepulsionDiffusionGradientReductionFilter<TInputScalarType, TOutputScalarType>
::GenerateData()
{
unsigned int randSeed = time(NULL);
if(m_InputBValueMap.empty() || m_NumGradientDirections.size()!=m_InputBValueMap.size())
return;
BValueMap manipulatedMap = m_InputBValueMap;
int shellCounter = 0;
for(BValueMap::iterator it = m_InputBValueMap.begin(); it != m_InputBValueMap.end(); it++ )
{
srand(randSeed);
// initialize index vectors
m_UsedGradientIndices.clear();
m_UnusedGradientIndices.clear();
if ( it->second.size() <= m_NumGradientDirections[shellCounter] )
{
itkWarningMacro( << "current directions: " << it->second.size() << " wanted directions: " << m_NumGradientDirections[shellCounter]);
m_NumGradientDirections[shellCounter] = it->second.size();
shellCounter++;
continue;
}
MITK_INFO << "Shell number: " << shellCounter;
unsigned int c=0;
for (unsigned int i=0; i<it->second.size(); i++)
{
if (c<m_NumGradientDirections[shellCounter])
m_UsedGradientIndices.push_back(it->second.at(i));
else
m_UnusedGradientIndices.push_back(it->second.at(i));
c++;
}
double minAngle = Costs();
double newMinAngle = 0;
MITK_INFO << "minimum angle: " << 180*minAngle/M_PI;
int stagnationCount = 0;
int rejectionCount = 0;
int maxRejections = m_NumGradientDirections[shellCounter] * 1000;
if (maxRejections<10000)
maxRejections = 10000;
int iUsed = 0;
if (m_UsedGradientIndices.size()>0)
while ( stagnationCount<1000 && rejectionCount<maxRejections )
{
// make proposal for new gradient configuration by randomly removing one of the currently used directions and instead adding one of the unused directions
//int iUsed = rand() % m_UsedGradientIndices.size();
int iUnUsed = rand() % m_UnusedGradientIndices.size();
int vUsed = m_UsedGradientIndices.at(iUsed);
int vUnUsed = m_UnusedGradientIndices.at(iUnUsed);
m_UsedGradientIndices.at(iUsed) = vUnUsed;
m_UnusedGradientIndices.at(iUnUsed) = vUsed;
newMinAngle = Costs(); // calculate costs of proposed configuration
if (newMinAngle > minAngle) // accept or reject proposal
{
MITK_INFO << "minimum angle: " << 180*newMinAngle/M_PI;
if ( (newMinAngle-minAngle)<0.01 )
stagnationCount++;
else
stagnationCount = 0;
minAngle = newMinAngle;
rejectionCount = 0;
}
else
{
rejectionCount++;
m_UsedGradientIndices.at(iUsed) = vUsed;
m_UnusedGradientIndices.at(iUnUsed) = vUnUsed;
}
iUsed++;
iUsed = iUsed % m_UsedGradientIndices.size();
}
manipulatedMap[it->first] = m_UsedGradientIndices;
shellCounter++;
}
int vecLength = 0 ;
for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++)
vecLength += it->second.size();
// initialize output image
typename OutputImageType::Pointer outImage = OutputImageType::New();
outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion());
outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->SetVectorLength( vecLength ); // Set the vector length
outImage->Allocate();
itk::ImageRegionIterator< OutputImageType > newIt(outImage, outImage->GetLargestPossibleRegion());
newIt.GoToBegin();
typename InputImageType::Pointer inImage = const_cast<InputImageType*>(this->GetInput(0));
itk::ImageRegionIterator< InputImageType > oldIt(inImage, inImage->GetLargestPossibleRegion());
oldIt.GoToBegin();
// initial new value of voxel
OutputPixelType newVec;
newVec.SetSize( vecLength );
newVec.AllocateElements( vecLength );
// generate new pixel values
while(!newIt.IsAtEnd())
{
// init new vector with zeros
newVec.Fill(0.0);
InputPixelType oldVec = oldIt.Get();
int index = 0;
for(BValueMap::iterator it=manipulatedMap.begin(); it!=manipulatedMap.end(); it++)
for(unsigned int j=0; j<it->second.size(); j++)
{
newVec[index] = oldVec[it->second.at(j)];
index++;
}
newIt.Set(newVec);
++newIt;
++oldIt;
}
// set new gradient directions
m_GradientDirections = GradientDirectionContainerType::New();
int index = 0;
for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++)
for(unsigned int j = 0; j < it->second.size(); j++)
{
m_GradientDirections->InsertElement(index, m_OriginalGradientDirections->at(it->second.at(j)));
index++;
}
this->SetNumberOfRequiredOutputs (1);
this->SetNthOutput (0, outImage);
MITK_INFO << "...done";
}
+template <class TInputScalarType, class TOutputScalarType>
+void
+ElectrostaticRepulsionDiffusionGradientReductionFilter<TInputScalarType, TOutputScalarType>
+::UpdateOutputInformation()
+{
+ Superclass::UpdateOutputInformation();
+ int vecLength = 0 ;
+ for(unsigned int index = 0; index < m_NumGradientDirections.size(); index++)
+ {
+ vecLength += m_NumGradientDirections[index];
+ }
+
+ this->GetOutput()->SetVectorLength( vecLength );
+}
} // end of namespace
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h
index 5db210b559..54fcc4e434 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h
@@ -1,78 +1,77 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
/*===================================================================
This file is based heavily on a corresponding ITK filter.
===================================================================*/
#ifndef __itkExtractDwiChannelFilter_h_
#define __itkExtractDwiChannelFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
-#include <mitkDiffusionImage.h>
#include <itkImageRegionIteratorWithIndex.h>
-#include <mitkDiffusionImage.h>
+#include <mitkImage.h>
namespace itk{
/** \class ExtractDwiChannelFilter
- * \brief Remove spcified channels from diffusion-weighted image.
+ * \brief Remove specified channels from diffusion-weighted image.
*/
template< class TInPixelType >
class ExtractDwiChannelFilter :
public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > >
{
public:
typedef ExtractDwiChannelFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > > Superclass;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(ExtractDwiChannelFilter, ImageToImageFilter)
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
itkSetMacro( ChannelIndex, unsigned int )
protected:
ExtractDwiChannelFilter();
~ExtractDwiChannelFilter() {}
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType );
unsigned int m_ChannelIndex;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkExtractDwiChannelFilter.txx"
#endif
#endif //__itkExtractDwiChannelFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h
index 878443938b..34cf4ef219 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h
@@ -1,147 +1,145 @@
/*===================================================================
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 __itkNonLocalMeansDenoisingFilter_h_
#define __itkNonLocalMeansDenoisingFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
-#include <mitkDiffusionImage.h>
-
namespace itk{
/** @class NonLocalMeansDenoisingFilter
* @brief This class denoises a vectorimage according to the non-local means procedure.
*
* This Filter needs as an input a diffusion weigthed image, which will be denoised unsing the non-local means principle.
* An input mask is optional to denoise only inside the mask range. All other voxels will be set to 0.
*/
template< class TPixelType >
class NonLocalMeansDenoisingFilter :
public ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > >
{
public:
/** Typedefs */
typedef NonLocalMeansDenoisingFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > > Superclass;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
typedef Image <TPixelType, 3> MaskImageType;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(NonLocalMeansDenoisingFilter, ImageToImageFilter)
/**
* @brief Set flag to use joint information
*/
itkSetMacro(UseJointInformation, bool)
/**
* @brief Set the searchradius
*
* The searchradius generates a neighborhood of size (2 * searchradius + 1)³.
* Default is 4.
*/
itkSetMacro(SearchRadius, int)
/**
* @brief Set the comparisonradius
*
* The comparisonradius generates neighborhoods of size (2 * comparisonradius +1)³.
* Default is 1.
*/
itkSetMacro(ComparisonRadius, int)
/**
* @brief Set the variance of the noise
*
* The variance of the noise needs to be estimated to use this filter properly.
* Default is 1.
*/
itkSetMacro(Variance, double)
/**
* @brief Set flag to use a rician adaption
*
* If this flag is true the filter uses a method which is optimized for Rician distributed noise.
*/
itkSetMacro(UseRicianAdaption, bool)
/**
* @brief Get the amount of calculated Voxels
*
* @return the number of calculated Voxels until yet, useful for the use of a progressbars.
*/
itkGetMacro(CurrentVoxelCount, unsigned int)
/** @brief Set the input image. **/
void SetInputImage(const InputImageType* image);
/**
* @brief Set a denoising mask
*
* optional
*
* Set a mask to denoise only the masked area, all voxel outside this area will be set to 0.
*/
void SetInputMask(MaskImageType* mask);
protected:
NonLocalMeansDenoisingFilter();
~NonLocalMeansDenoisingFilter() {}
/**
* @brief Calculations which need to be done before the denoising starts
*
* This method is called before the denoising starts. It calculates the ROI if a mask is used
* and sets the number of processed voxels to zero.
*/
void BeforeThreadedGenerateData();
/**
* @brief Denoising procedure
*
* This method calculates the denoised voxelvalue for each voxel in the image in multiple threads.
* If a mask is used, voxels outside the masked area will be set to 0.
*
* @param outputRegionForThread Region to denoise for each thread.
*/
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType);
private:
int m_SearchRadius; ///< Radius of the searchblock.
int m_ComparisonRadius; ///< Radius of the comparisonblock.
bool m_UseJointInformation; ///< Flag to use joint information.
bool m_UseRicianAdaption; ///< Flag to use rician adaption.
unsigned int m_CurrentVoxelCount; ///< Amount of processed voxels.
double m_Variance; ///< Estimated noise variance.
typename MaskImageType::Pointer m_Mask; ///< Pointer to the mask image.
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkNonLocalMeansDenoisingFilter.txx"
#endif
#endif //__itkNonLocalMeansDenoisingFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h
index 7e18435d36..168a92b5f8 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h
@@ -1,85 +1,84 @@
/*===================================================================
The Medical Imaging Interaction Toolkit (MITK)
Copyright (c) German Cancer Research Center,
Division of Medical and Biological Informatics.
All rights reserved.
This software is distributed WITHOUT ANY WARRANTY; without
even the implied warranty of MERCHANTABILITY or FITNESS FOR
A PARTICULAR PURPOSE.
See LICENSE.txt or http://www.mitk.org for details.
===================================================================*/
/*===================================================================
This file is based heavily on a corresponding ITK filter.
===================================================================*/
#ifndef __itkRemoveDwiChannelFilter_h_
#define __itkRemoveDwiChannelFilter_h_
#include "itkImageToImageFilter.h"
#include "itkVectorImage.h"
-#include <mitkDiffusionImage.h>
#include <itkImageRegionIteratorWithIndex.h>
-#include <mitkDiffusionImage.h>
+#include <mitkDiffusionPropertyHelper.h>
namespace itk{
/** \class RemoveDwiChannelFilter
* \brief Remove spcified channels from diffusion-weighted image.
*/
template< class TInPixelType >
class RemoveDwiChannelFilter :
public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > >
{
public:
typedef RemoveDwiChannelFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(RemoveDwiChannelFilter, ImageToImageFilter)
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename Superclass::OutputImageRegionType OutputImageRegionType;
- typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionType DirectionType;
- typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionContainerType DirectionContainerType;
+ typedef typename mitk::DiffusionPropertyHelper::GradientDirectionType DirectionType;
+ typedef typename mitk::DiffusionPropertyHelper::GradientDirectionsContainerType DirectionContainerType;
void SetChannelIndices( std::vector< unsigned int > indices ){ m_ChannelIndices = indices; }
void SetDirections( typename DirectionContainerType::Pointer directions ){ m_Directions = directions; }
typename DirectionContainerType::Pointer GetNewDirections(){ return m_NewDirections; }
protected:
RemoveDwiChannelFilter();
~RemoveDwiChannelFilter() {}
void PrintSelf(std::ostream& os, Indent indent) const;
void BeforeThreadedGenerateData();
void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType id );
std::vector< unsigned int > m_ChannelIndices;
typename DirectionContainerType::Pointer m_Directions;
typename DirectionContainerType::Pointer m_NewDirections;
};
}
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkRemoveDwiChannelFilter.txx"
#endif
#endif //__itkRemoveDwiChannelFilter_h_
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h
index 2bbaf13288..925d1ff77a 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h
@@ -1,138 +1,146 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkResampleDwiImageFilter.h $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_ResampleDwiImageFilter_h_
#define _itk_ResampleDwiImageFilter_h_
#include <itkImageToImageFilter.h>
#include <itkVectorImage.h>
#include <itkPointShell.h>
namespace itk
{
/**
* \brief Resample DWI channel by channel. */
template <class TScalarType>
class ResampleDwiImageFilter
: public ImageToImageFilter<itk::VectorImage<TScalarType,3>, itk::VectorImage<TScalarType,3> >
{
public:
enum Interpolation {
Interpolate_NearestNeighbour,
Interpolate_Linear,
Interpolate_BSpline,
Interpolate_WindowedSinc
};
typedef ResampleDwiImageFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< itk::VectorImage<TScalarType,3>, itk::VectorImage<TScalarType,3> >
Superclass;
typedef itk::Vector< double, 3 > DoubleVectorType;
typedef itk::VectorImage<TScalarType,3> DwiImageType;
typedef itk::Image<TScalarType,3> DwiChannelType;
/** Method for creation through the object factory. */
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
/** Runtime information support. */
itkTypeMacro(ResampleDwiImageFilter, ImageToImageFilter)
itkSetMacro( Interpolation, Interpolation )
void SetSamplingFactor(DoubleVectorType sampling)
{
m_NewSpacing = this->GetInput()->GetSpacing();
m_NewSpacing[0] /= sampling[0];
m_NewSpacing[1] /= sampling[1];
m_NewSpacing[2] /= sampling[2];
m_NewImageRegion = this->GetInput()->GetLargestPossibleRegion();
m_NewImageRegion.SetSize(0, m_NewImageRegion.GetSize(0)*sampling[0]);
m_NewImageRegion.SetSize(1, m_NewImageRegion.GetSize(1)*sampling[1]);
m_NewImageRegion.SetSize(2, m_NewImageRegion.GetSize(2)*sampling[2]);
}
void SetNewSpacing(DoubleVectorType spacing)
{
DoubleVectorType oldSpacing = this->GetInput()->GetSpacing();
DoubleVectorType sampling;
sampling[0] = oldSpacing[0]/spacing[0];
sampling[1] = oldSpacing[1]/spacing[1];
sampling[2] = oldSpacing[2]/spacing[2];
m_NewSpacing = spacing;
m_NewImageRegion = this->GetInput()->GetLargestPossibleRegion();
m_NewImageRegion.SetSize(0, m_NewImageRegion.GetSize(0)*sampling[0]);
m_NewImageRegion.SetSize(1, m_NewImageRegion.GetSize(1)*sampling[1]);
m_NewImageRegion.SetSize(2, m_NewImageRegion.GetSize(2)*sampling[2]);
}
void SetNewImageSize(ImageRegion<3> region)
{
ImageRegion<3> oldRegion = this->GetInput()->GetLargestPossibleRegion();
DoubleVectorType sampling;
sampling[0] = (double)region.GetSize(0)/oldRegion.GetSize(0);
sampling[1] = (double)region.GetSize(1)/oldRegion.GetSize(1);
sampling[2] = (double)region.GetSize(2)/oldRegion.GetSize(2);
m_NewImageRegion = region;
m_NewSpacing = this->GetInput()->GetSpacing();
m_NewSpacing[0] /= sampling[0];
m_NewSpacing[1] /= sampling[1];
m_NewSpacing[2] /= sampling[2];
}
+ virtual void UpdateOutputInformation();
+
+ virtual void PropagateRequestedRegion(){}
+
+ virtual void PropagateRequestedRegion(itk::DataObject *output){}
+
+ virtual void VerifyInputInformation(){}
+
protected:
ResampleDwiImageFilter();
~ResampleDwiImageFilter(){}
void GenerateData();
DoubleVectorType m_NewSpacing;
ImageRegion<3> m_NewImageRegion;
Interpolation m_Interpolation;
};
} // end of namespace
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkResampleDwiImageFilter.txx"
#endif
#endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx
index 516c121104..b317f88f4d 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx
@@ -1,166 +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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkResampleDwiImageFilter.txx $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_ResampleDwiImageFilter_txx_
#define _itk_ResampleDwiImageFilter_txx_
#endif
#define _USE_MATH_DEFINES
#include "itkResampleDwiImageFilter.h"
#include <itkImageRegionIterator.h>
#include <itkImageRegion.h>
#include <itkResampleImageFilter.h>
#include <itkBSplineInterpolateImageFunction.h>
#include <itkNearestNeighborInterpolateImageFunction.h>
#include <itkLinearInterpolateImageFunction.h>
#include <itkWindowedSincInterpolateImageFunction.h>
namespace itk
{
template <class TScalarType>
ResampleDwiImageFilter<TScalarType>
::ResampleDwiImageFilter()
: m_Interpolation(Interpolate_Linear)
{
this->SetNumberOfRequiredInputs( 1 );
}
template <class TScalarType>
void
ResampleDwiImageFilter<TScalarType>
::GenerateData()
{
// // initialize output image
// itk::Vector< double, 3 > spacing = this->GetInput()->GetSpacing();
// spacing[0] /= m_SamplingFactor[0];
// spacing[1] /= m_SamplingFactor[1];
// spacing[2] /= m_SamplingFactor[2];
// ImageRegion<3> region = this->GetInput()->GetLargestPossibleRegion();
// region.SetSize(0, region.GetSize(0)*m_SamplingFactor[0]);
// region.SetSize(1, region.GetSize(1)*m_SamplingFactor[1]);
// region.SetSize(2, region.GetSize(2)*m_SamplingFactor[2]);
itk::Point<double,3> origin = this->GetInput()->GetOrigin();
origin[0] -= this->GetInput()->GetSpacing()[0]/2;
origin[1] -= this->GetInput()->GetSpacing()[1]/2;
origin[2] -= this->GetInput()->GetSpacing()[2]/2;
origin[0] += m_NewSpacing[0]/2;
origin[1] += m_NewSpacing[1]/2;
origin[2] += m_NewSpacing[2]/2;
typename DwiImageType::Pointer outImage = DwiImageType::New();
outImage->SetSpacing( m_NewSpacing );
outImage->SetOrigin( origin );
outImage->SetDirection( this->GetInput()->GetDirection() );
outImage->SetLargestPossibleRegion( m_NewImageRegion );
outImage->SetBufferedRegion( m_NewImageRegion );
outImage->SetRequestedRegion( m_NewImageRegion );
outImage->SetVectorLength( this->GetInput()->GetVectorLength() );
outImage->Allocate();
typename itk::ResampleImageFilter<DwiChannelType, DwiChannelType>::Pointer resampler = itk::ResampleImageFilter<DwiChannelType, DwiChannelType>::New();
resampler->SetOutputParametersFromImage(outImage);
switch (m_Interpolation)
{
case Interpolate_NearestNeighbour:
{
typename itk::NearestNeighborInterpolateImageFunction<DwiChannelType>::Pointer interp = itk::NearestNeighborInterpolateImageFunction<DwiChannelType>::New();
resampler->SetInterpolator(interp);
break;
}
case Interpolate_Linear:
{
typename itk::LinearInterpolateImageFunction<DwiChannelType>::Pointer interp = itk::LinearInterpolateImageFunction<DwiChannelType>::New();
resampler->SetInterpolator(interp);
break;
}
case Interpolate_BSpline:
{
typename itk::BSplineInterpolateImageFunction<DwiChannelType>::Pointer interp = itk::BSplineInterpolateImageFunction<DwiChannelType>::New();
resampler->SetInterpolator(interp);
break;
}
case Interpolate_WindowedSinc:
{
typename itk::WindowedSincInterpolateImageFunction<DwiChannelType, 3>::Pointer interp = itk::WindowedSincInterpolateImageFunction<DwiChannelType, 3>::New();
resampler->SetInterpolator(interp);
break;
}
default:
{
typename itk::LinearInterpolateImageFunction<DwiChannelType>::Pointer interp = itk::LinearInterpolateImageFunction<DwiChannelType>::New();
resampler->SetInterpolator(interp);
}
}
for (unsigned int i=0; i<this->GetInput()->GetVectorLength(); i++)
{
typename DwiChannelType::Pointer channel = DwiChannelType::New();
channel->SetSpacing( this->GetInput()->GetSpacing() );
channel->SetOrigin( this->GetInput()->GetOrigin() );
channel->SetDirection( this->GetInput()->GetDirection() );
channel->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion() );
channel->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
channel->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
channel->Allocate();
ImageRegionIterator<DwiChannelType> it(channel, channel->GetLargestPossibleRegion());
while(!it.IsAtEnd())
{
typename DwiImageType::PixelType pix = this->GetInput()->GetPixel(it.GetIndex());
it.Set(pix[i]);
++it;
}
resampler->SetInput(channel);
resampler->Update();
channel = resampler->GetOutput();
ImageRegionIterator<DwiImageType> it2(outImage, outImage->GetLargestPossibleRegion());
while(!it2.IsAtEnd())
{
typename DwiImageType::PixelType pix = it2.Get();
pix[i] = channel->GetPixel(it2.GetIndex());
it2.Set(pix);
++it2;
}
}
this->SetNthOutput(0, outImage);
}
+template <class TScalarType>
+void
+ResampleDwiImageFilter<TScalarType>
+::UpdateOutputInformation()
+{
+ // Calls to superclass updateoutputinformation
+ //Superclass::UpdateOutputInformation();
+ this->GetOutput()->SetSpacing(m_NewSpacing);
+ this->GetOutput()->SetLargestPossibleRegion(m_NewImageRegion);
+}
} // end of namespace
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx
index 49715a1226..fb1c06a99a 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx
@@ -1,285 +1,285 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_ResidualImageFilter_txx_
#define _itk_ResidualImageFilter_txx_
#endif
#include "itkResidualImageFilter.h"
#include <mitkCommon.h>
#include <cmath>
namespace itk
{
template <class TInputScalarType, class TOutputScalarType>
void
ResidualImageFilter<TInputScalarType, TOutputScalarType>
::GenerateData()
{
typename InputImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize();
typename InputImageType::SizeType size2 = m_SecondDiffusionImage->GetLargestPossibleRegion().GetSize();
if(size != size2)
{
MITK_ERROR << "Sizes do not match";
return;
}
// Initialize output image
typename OutputImageType::Pointer outputImage =
static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
outputImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
outputImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
outputImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
outputImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() );
outputImage->Allocate();
outputImage->FillBuffer(0.0);
std::vector< std::vector<double> > residuals;
// per slice, per volume
std::vector< std::vector <std::vector<double> > > residualsPerSlice;
// Detrmine number of B0 images
int numberB0=0;
for(unsigned int i=0; i<m_Gradients->Size(); i++)
{
GradientDirectionType grad = m_Gradients->ElementAt(i);
if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001)
{
numberB0++;
}
}
residuals.resize(this->GetInput()->GetVectorLength()-numberB0);
// Calculate the standard residual image and for each volume put all residuals in a vector
for(unsigned int z=0; z<size[2]; z++)
{
std::vector< std::vector<double> > sliceResiduals; // residuals per volume for this slice
sliceResiduals.resize(this->GetInput()->GetVectorLength()-numberB0);
for(unsigned int y=0; y<size[1]; y++)
{
for(unsigned int x=0; x<size[0]; x++)
{
// Check if b0 exceeds threshold
itk::Index<3> ix;
ix[0] = x;
ix[1] = y;
ix[2] = z;
typename InputImageType::PixelType p1 = this->GetInput()->GetPixel(ix);
typename InputImageType::PixelType p2 = m_SecondDiffusionImage->GetPixel(ix);
if(p1.GetSize() != p2.GetSize())
{
MITK_ERROR << "Vector sizes do not match";
return;
}
if(p1.GetElement(m_B0Index) <= m_B0Threshold)
{
continue;
}
double res = 0;
int shift = 0; // correction for the skipped B0 images
for(unsigned int i = 0; i<p1.GetSize(); i++)
{
GradientDirectionType grad = m_Gradients->ElementAt(i);
if(!(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001))
{
double val1 = (double)p1.GetElement(i);
double val2 = (double)p2.GetElement(i);
res += abs(val1-val2);
residuals[i-shift].push_back(val1-val2);
sliceResiduals[i-shift].push_back(val1-val2);
}
else
{
shift++;
}
}
res = res/p1.GetSize();
outputImage->SetPixel(ix, res);
}
}
residualsPerSlice.push_back(sliceResiduals);
}
// for each dw volume: sort the the measured residuals (for each voxel) to enable determining Q1 and Q3; calculate means
// determine percentage of errors as described in QUALITY ASSESSMENT THROUGH ANALYSIS OF RESIDUALS OF DIFFUSION IMAGE FITTING
// Leemans et al 2008
double q1,q3, median;
std::vector< std::vector<double> >::iterator it = residuals.begin();
while(it != residuals.end())
{
std::vector<double> res = *it;
// sort
std::sort(res.begin(), res.end());
q1 = res[0.25*res.size()];
m_Q1.push_back(q1);
q3 = res[0.75*res.size()];
m_Q3.push_back(q3);
median = res[0.5*res.size()];
double iqr = q3-q1;
double outlierThreshold = median + 1.5*iqr;
double numberOfOutliers = 0.0;
std::vector<double>::iterator resIt = res.begin();
double mean = 0;
while(resIt != res.end())
{
double f = *resIt;
if(f>outlierThreshold)
{
numberOfOutliers++;
}
mean += f;
++resIt;
}
double percOfOutliers = 100 * numberOfOutliers / res.size();
m_PercentagesOfOutliers.push_back(percOfOutliers);
mean /= res.size();
m_Means.push_back(mean);
++it;
}
// Calculate for each slice the number of outliers per volume(dw volume)
std::vector< std::vector <std::vector<double> > >::iterator sliceIt = residualsPerSlice.begin();
while(sliceIt != residualsPerSlice.end())
{
std::vector< std::vector<double> > currentSlice = *sliceIt;
std::vector<double> percentages;
std::vector< std::vector<double> >::iterator volIt = currentSlice.begin();
while(volIt != currentSlice.end())
{
std::vector<double> currentVolume = *volIt;
//sort
std::sort(currentVolume.begin(), currentVolume.end());
q1 = currentVolume[0.25*currentVolume.size()];
q3 = currentVolume[0.75*currentVolume.size()];
median = currentVolume[0.5*currentVolume.size()];
double iqr = q3-q1;
double outlierThreshold = median + 1.5*iqr;
double numberOfOutliers = 0.0;
std::vector<double>::iterator resIt = currentVolume.begin();
- double mean;
+ double mean(0.0);
while(resIt != currentVolume.end())
{
double f = *resIt;
if(f>outlierThreshold)
{
numberOfOutliers++;
}
mean += f;
++resIt;
}
double percOfOutliers = 100 * numberOfOutliers / currentVolume.size();
percentages.push_back(percOfOutliers);
++volIt;
}
m_OutliersPerSlice.push_back(percentages);
++sliceIt;
}
}
} // end of namespace
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx
index 34cb1f08cb..ad119b1a73 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx
@@ -1,145 +1,145 @@
/*===================================================================
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 __itkTensorDerivedMeasurementsFilter_txx
#define __itkTensorDerivedMeasurementsFilter_txx
-
+#include <itkImageRegionIterator.h>
namespace itk {
template <class TPixel>
TensorDerivedMeasurementsFilter<TPixel>::TensorDerivedMeasurementsFilter() : m_Measure(AD)
{
}
template <class TPixel>
void TensorDerivedMeasurementsFilter<TPixel>::GenerateData()
{
typename TensorImageType::Pointer tensorImage = static_cast< TensorImageType * >( this->ProcessObject::GetInput(0) );
typedef ImageRegionConstIterator< TensorImageType > TensorImageIteratorType;
typedef ImageRegionIterator< OutputImageType > OutputImageIteratorType;
typename OutputImageType::Pointer outputImage =
static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
typename TensorImageType::RegionType region = tensorImage->GetLargestPossibleRegion();
outputImage->SetSpacing( tensorImage->GetSpacing() ); // Set the image spacing
outputImage->SetOrigin( tensorImage->GetOrigin() ); // Set the image origin
outputImage->SetDirection( tensorImage->GetDirection() ); // Set the image direction
outputImage->SetRegions( tensorImage->GetLargestPossibleRegion());
outputImage->Allocate();
TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetLargestPossibleRegion());
OutputImageIteratorType outputIt(outputImage, outputImage->GetLargestPossibleRegion());
tensorIt.GoToBegin();
outputIt.GoToBegin();
while(!tensorIt.IsAtEnd() && !outputIt.IsAtEnd()){
TensorType tensor = tensorIt.Get();
switch(m_Measure)
{
case FA:
{
TPixel diffusionIndex = tensor.GetFractionalAnisotropy();
outputIt.Set(diffusionIndex);
break;
}
case RA:
{
TPixel diffusionIndex = tensor.GetRelativeAnisotropy();
outputIt.Set(diffusionIndex);
break;
}
case AD:
{
// eigenvalues are sorted in ascending order by default because the
// itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
outputIt.Set(evs[2]);
break;
}
case RD:
{
// eigenvalues are sorted in ascending order by default because the
// itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
outputIt.Set((evs[0]+evs[1])/2.0);
break;
}
case CA:
{
// eigenvalues are sorted in ascending order by default because the
// itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
if (evs[2] == 0)
{
outputIt.Set(0);
break;
}
outputIt.Set(1.0-(evs[0]+evs[1])/(2.0*evs[2]));
break;
}
case L2:
{
// eigenvalues are sorted in ascending order by default because the
// itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
outputIt.Set(evs[1]);
break;
}
case L3:
{
// eigenvalues are sorted in ascending order by default because the
// itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
outputIt.Set(evs[0]);
break;
}
case MD:
{
typename TensorType::EigenValuesArrayType evs;
tensor.ComputeEigenValues(evs);
outputIt.Set((evs[0]+evs[1]+evs[2])/3.0);
break;
}
}
++tensorIt;
++outputIt;
}
}
}
#endif // __itkTensorDerivedMeasurements_txx
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h
index cbed2d792b..b588be572a 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h
@@ -1,168 +1,171 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.h $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_TensorImageToDiffusionImageFilter_h_
#define _itk_TensorImageToDiffusionImageFilter_h_
#include "itkImageToImageFilter.h"
#include <itkDiffusionTensor3D.h>
#include <itkVectorContainer.h>
namespace itk
{
template <class TInputScalarType, class TOutputScalarType>
class TensorImageToDiffusionImageFilter
: public ImageToImageFilter<itk::Image<itk::DiffusionTensor3D<TInputScalarType>,3>, itk::VectorImage<TOutputScalarType,3> >
{
public:
typedef TInputScalarType InputScalarType;
typedef itk::DiffusionTensor3D<InputScalarType> InputPixelType;
typedef itk::Image<InputPixelType,3> InputImageType;
typedef typename InputImageType::RegionType InputImageRegionType;
typedef TOutputScalarType OutputScalarType;
typedef itk::VectorImage<OutputScalarType,3> OutputImageType;
typedef typename OutputImageType::PixelType OutputPixelType;
typedef typename OutputImageType::RegionType OutputImageRegionType;
typedef OutputScalarType BaselineScalarType;
typedef BaselineScalarType BaselinePixelType;
typedef typename itk::Image<BaselinePixelType,3> BaselineImageType;
typedef typename BaselineImageType::RegionType BaselineImageRegionType;
typedef itk::Image< short, 3> MaskImageType;
typedef MaskImageType::RegionType MaskImageRegionType;
typedef TensorImageToDiffusionImageFilter Self;
typedef ImageToImageFilter<InputImageType, OutputImageType> Superclass;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
itkTypeMacro (TensorImageToDiffusionImageFilter, ImageToImageFilter);
itkStaticConstMacro (ImageDimension, unsigned int,
OutputImageType::ImageDimension);
itkFactorylessNewMacro(Self)
itkCloneMacro(Self)
typedef vnl_vector_fixed<double,3> GradientType;
typedef VectorContainer<unsigned int, GradientType> GradientListType;
typedef GradientListType::Pointer GradientListPointerType;
/** Manually Set/Get a list of gradients */
void SetGradientList(const GradientListPointerType list)
{
m_GradientList = list;
this->Modified();
}
GradientListPointerType GetGradientList(void) const
{return m_GradientList;}
void SetBValue( const double& bval)
{ m_BValue = bval; }
void SetMaskImage( MaskImageType::Pointer maskimage )
{
m_MaskImage = maskimage;
}
/**
* @brief Set an external baseline image for signal generation (optional)
*
* An option to enforce a specific baseline image. If none provided (default) the filter uses
* the itk::TensorToL2NormImageFilter to generate the modelled baseline image.
*/
void SetExternalBaselineImage( typename BaselineImageType::Pointer bimage)
{
m_BaselineImage = bimage;
}
itkSetMacro(Min, OutputScalarType);
itkSetMacro(Max, OutputScalarType);
protected:
TensorImageToDiffusionImageFilter()
{
m_BValue = 1.0;
m_BaselineImage = 0;
m_Min = 0.0;
m_Max = 10000.0;
}
virtual ~TensorImageToDiffusionImageFilter(){}
void PrintSelf (std::ostream& os, Indent indent) const
{
Superclass::PrintSelf (os, indent);
}
virtual void BeforeThreadedGenerateData( void );
virtual void ThreadedGenerateData( const
OutputImageRegionType &outputRegionForThread, ThreadIdType);
//void GenerateData();
- GradientListPointerType m_GradientList;
+ virtual void UpdateOutputInformation();
+
+
+ private:
+
+ TensorImageToDiffusionImageFilter (const Self&);
+ void operator=(const Self&);
+
+ GradientListType::Pointer m_GradientList;
double m_BValue;
typename BaselineImageType::Pointer m_BaselineImage;
OutputScalarType m_Min;
OutputScalarType m_Max;
MaskImageType::Pointer m_MaskImage;
- private:
-
- TensorImageToDiffusionImageFilter (const Self&);
- void operator=(const Self&);
-
};
} // end of namespace
#ifndef ITK_MANUAL_INSTANTIATION
#include "itkTensorImageToDiffusionImageFilter.txx"
#endif
#endif
diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx
index e812dfb51a..67a9b355d2 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx
+++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx
@@ -1,235 +1,244 @@
/*===================================================================
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.
===================================================================*/
/*=========================================================================
2
3 Program: Tensor ToolKit - TTK
4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $
5 Language: C++
6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
7 Version: $Revision: 68 $
8
9 Copyright (c) INRIA 2010. All rights reserved.
10 See LICENSE.txt for details.
11
12 This software is distributed WITHOUT ANY WARRANTY; without even
13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14 PURPOSE. See the above copyright notices for more information.
15
16 =========================================================================*/
#ifndef _itk_TensorImageToDiffusionImageFilter_txx_
#define _itk_TensorImageToDiffusionImageFilter_txx_
#endif
#include "itkTensorImageToDiffusionImageFilter.h"
#include "itkTensorToL2NormImageFilter.h"
#include "itkRescaleIntensityImageFilter.h"
#include <itkImageRegionIterator.h>
#include <itkImageRegionConstIterator.h>
namespace itk
{
template <class TInputScalarType, class TOutputScalarType>
void
TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
::BeforeThreadedGenerateData()
{
if( m_GradientList->Size()==0 )
{
throw itk::ExceptionObject (__FILE__,__LINE__,"Error: gradient list is empty, cannot generate DWI.");
}
if( m_BaselineImage.IsNull() )
{
// create a B0 image by taking the norm of the tensor field * scale:
typedef itk::TensorToL2NormImageFilter<InputImageType, itk::Image<InputScalarType,3> >
TensorToL2NormFilterType;
typename TensorToL2NormFilterType::Pointer myFilter1 = TensorToL2NormFilterType::New();
myFilter1->SetInput (this->GetInput());
try
{
myFilter1->Update();
}
catch (itk::ExceptionObject &e)
{
std::cerr << e;
return;
}
typename itk::RescaleIntensityImageFilter< itk::Image<InputScalarType,3>, BaselineImageType>::Pointer rescaler=
itk::RescaleIntensityImageFilter<itk::Image<InputScalarType,3>, BaselineImageType>::New();
rescaler->SetOutputMinimum ( m_Min );
rescaler->SetOutputMaximum ( m_Max );
rescaler->SetInput ( myFilter1->GetOutput() );
try
{
rescaler->Update();
}
catch (itk::ExceptionObject &e)
{
std::cerr << e;
return;
}
m_BaselineImage = rescaler->GetOutput();
}
typename OutputImageType::Pointer outImage = OutputImageType::New();
outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion());
outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->SetVectorLength(m_GradientList->Size());
outImage->Allocate();
this->SetNumberOfRequiredOutputs (1);
this->SetNthOutput (0, outImage);
}
template <class TInputScalarType, class TOutputScalarType>
void TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId )
{
typedef ImageRegionIterator<OutputImageType> IteratorOutputType;
typedef ImageRegionConstIterator<InputImageType> IteratorInputType;
typedef ImageRegionConstIterator<BaselineImageType> IteratorBaselineType;
unsigned long numPixels = outputRegionForThread.GetNumberOfPixels();
unsigned long step = numPixels/100;
unsigned long progress = 0;
IteratorOutputType itOut (this->GetOutput(), outputRegionForThread);
IteratorInputType itIn (this->GetInput(), outputRegionForThread);
IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread);
typedef ImageRegionConstIterator< MaskImageType > IteratorMaskImageType;
IteratorMaskImageType itMask;
if( m_MaskImage.IsNotNull() )
{
itMask = IteratorMaskImageType( m_MaskImage, outputRegionForThread);
itMask.GoToBegin();
}
if( threadId==0 )
{
this->UpdateProgress (0.0);
}
while(!itIn.IsAtEnd())
{
if( this->GetAbortGenerateData() )
{
throw itk::ProcessAborted(__FILE__,__LINE__);
}
InputPixelType T = itIn.Get();
BaselinePixelType b0 = itB0.Get();
OutputPixelType out;
out.SetSize(m_GradientList->Size());
out.Fill(0);
short maskvalue = 1;
if( m_MaskImage.IsNotNull() )
{
maskvalue = itMask.Get();
++itMask;
}
std::vector<unsigned int> b0_indices;
if( b0 > 0)
{
for( unsigned int i=0; i<m_GradientList->Size(); i++)
{
GradientType g = m_GradientList->at(i);
// normalize vector so the following computations work
const double twonorm = g.two_norm();
if( twonorm < vnl_math::eps )
{
b0_indices.push_back(i);
continue;
}
GradientType gn = g.normalize();
InputPixelType S;
S[0] = gn[0]*gn[0];
S[1] = gn[1]*gn[0];
S[2] = gn[2]*gn[0];
S[3] = gn[1]*gn[1];
S[4] = gn[2]*gn[1];
S[5] = gn[2]*gn[2];
const double res =
T[0]*S[0] +
2 * T[1]*S[1] + T[3]*S[3] +
2 * T[2]*S[2] + 2 * T[4]*S[4] + T[5]*S[5];
// check for corrupted tensor
if (res>=0)
{
// estimate the bvalue from the base value and the norm of the gradient
// - because of this estimation the vector have to be normalized beforehand
// otherwise the modelled signal is wrong ( i.e. not scaled properly )
const double bval = m_BValue * twonorm * twonorm;
out[i] = static_cast<OutputScalarType>( maskvalue * 1.0 * b0 * exp ( -1.0 * bval * res ) );
}
}
}
for(unsigned int idx = 0; idx < b0_indices.size(); idx++ )
{
out[b0_indices.at(idx)] = b0;
}
itOut.Set(out);
if( threadId==0 && step>0)
{
if( (progress%step)==0 )
{
this->UpdateProgress ( double(progress)/double(numPixels) );
}
}
++progress;
++itB0;
++itIn;
++itOut;
}
if( threadId==0 )
{
this->UpdateProgress (1.0);
}
}
+template <class TInputScalarType, class TOutputScalarType>
+void
+TensorImageToDiffusionImageFilter<TInputScalarType, TOutputScalarType>
+::UpdateOutputInformation()
+{
+ // Calls to superclass updateoutputinformation
+ Superclass::UpdateOutputInformation();
+ this->GetOutput()->SetVectorLength( m_GradientList->size() );
+}
} // end of namespace
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp
index b4fb13c110..93dc5802af 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp
@@ -1,306 +1,311 @@
/*===================================================================
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 "mitkDiffusionDICOMFileReader.h"
#include "mitkDiffusionDICOMFileReaderHelper.h"
#include "mitkDiffusionHeaderSiemensDICOMFileReader.h"
#include "mitkDiffusionHeaderSiemensMosaicDICOMFileReader.h"
#include "mitkDiffusionHeaderGEDICOMFileReader.h"
#include "mitkDiffusionHeaderPhilipsDICOMFileReader.h"
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkITKImageImport.h>
#include "mitkStringProperty.h"
+#include <mitkImageCast.h>
static void PerformHeaderAnalysis( mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType headers )
{
unsigned int images = headers.size();
unsigned int unweighted_images = 0;
unsigned int weighted_images = 0;
mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType::const_iterator c_iter = headers.begin();
while( c_iter != headers.end() )
{
const mitk::DiffusionImageDICOMHeaderInformation h = *c_iter;
if( h.baseline ) unweighted_images++;
if( h.b_value > 0 ) weighted_images++;
++c_iter;
}
MITK_INFO << " :: Analyzed volumes " << images << "\n"
<< " :: \t"<< unweighted_images << " b = 0" << "\n"
<< " :: \t"<< weighted_images << " b > 0";
}
mitk::DiffusionDICOMFileReader::DiffusionDICOMFileReader()
{
}
mitk::DiffusionDICOMFileReader::~DiffusionDICOMFileReader()
{
}
bool mitk::DiffusionDICOMFileReader
::LoadImages()
{
unsigned int numberOfOutputs = this->GetNumberOfOutputs();
bool success = true;
for(unsigned int o = 0; o < numberOfOutputs; ++o)
{
success &= this->LoadSingleOutputImage( this->m_OutputHeaderContainer.at(o),
this->InternalGetOutput(o), this->m_IsMosaicData.at(o) );
}
return success;
}
bool mitk::DiffusionDICOMFileReader
::LoadSingleOutputImage( DiffusionHeaderDICOMFileReader::DICOMHeaderListType retrievedHeader,
DICOMImageBlockDescriptor& block, bool is_mosaic)
{
// prepare data reading
DiffusionDICOMFileReaderHelper helper;
DiffusionDICOMFileReaderHelper::VolumeFileNamesContainer filenames;
const DICOMImageFrameList& frames = block.GetImageFrameList();
int numberOfDWImages = block.GetIntProperty("timesteps", 1);
int numberOfFramesPerDWImage = frames.size() / numberOfDWImages;
assert( int( double((double) frames.size() / (double) numberOfDWImages)) == numberOfFramesPerDWImage );
for( int idx = 0; idx < numberOfDWImages; idx++ )
{
std::vector< std::string > FileNamesPerVolume;
DICOMImageFrameList::const_iterator timeStepStart = frames.begin() + idx * numberOfFramesPerDWImage;
DICOMImageFrameList::const_iterator timeStepEnd = frames.begin() + (idx+1) * numberOfFramesPerDWImage;
for (DICOMImageFrameList::const_iterator frameIter = timeStepStart;
frameIter != timeStepEnd;
++frameIter)
{
FileNamesPerVolume.push_back( (*frameIter)->Filename );
}
filenames.push_back( FileNamesPerVolume );
}
// TODO : only prototyping to test loading of diffusion images
// we need some solution for the different types
- typedef mitk::DiffusionImage<short> DiffusionImageType;
- DiffusionImageType::Pointer output_image = DiffusionImageType::New();
+ mitk::Image::Pointer output_image = mitk::Image::New();
- DiffusionImageType::GradientDirectionContainerType::Pointer directions =
- DiffusionImageType::GradientDirectionContainerType::New();
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer directions =
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New();
double max_bvalue = 0;
for( int idx = 0; idx < numberOfDWImages; idx++ )
{
DiffusionImageDICOMHeaderInformation header = retrievedHeader.at(idx);
if( max_bvalue < header.b_value )
max_bvalue = header.b_value;
}
// normalize the retrieved gradient directions according to the set b-value (maximal one)
for( int idx = 0; idx < numberOfDWImages; idx++ )
{
DiffusionImageDICOMHeaderInformation header = retrievedHeader.at(idx);
- DiffusionImageType::GradientDirectionType grad = header.g_vector;
+ mitk::DiffusionPropertyHelper::GradientDirectionType grad = header.g_vector;
grad.normalize();
grad *= sqrt( header.b_value / max_bvalue );
directions->push_back( grad );
}
// initialize the output image
- output_image->SetReferenceBValue( max_bvalue );
- output_image->SetDirections( directions );
+ output_image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( directions ) );
+ output_image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( max_bvalue ) );
if( is_mosaic && this->m_ResolveMosaic )
{
mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::Pointer mosaic_reader =
mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::New();
// retrieve the remaining meta-information needed for mosaic reconstruction
// it suffices to get it exemplatory from the first file in the file list
mosaic_reader->RetrieveMosaicInformation( filenames.at(0).at(0) );
mitk::MosaicDescriptor mdesc = mosaic_reader->GetMosaicDescriptor();
- output_image->SetVectorImage( helper.LoadMosaicToVector<short, 3>( filenames, mdesc ) );
+ mitk::CastToMitkImage( helper.LoadMosaicToVector<short, 3>( filenames, mdesc ), output_image );
}
else
{
- output_image->SetVectorImage( helper.LoadToVector<short, 3>( filenames ) );
+ mitk::CastToMitkImage( helper.LoadToVector<short, 3>( filenames ), output_image );
}
- output_image->InitializeFromVectorImage();
+
+ mitk::DiffusionPropertyHelper propertyHelper( output_image );
+ propertyHelper.InitializeImage();
+
output_image->SetProperty("diffusion.dicom.importname", mitk::StringProperty::New( helper.GetOutputName(filenames) ) );
block.SetMitkImage( (mitk::Image::Pointer) output_image );
return block.GetMitkImage().IsNotNull();
}
void mitk::DiffusionDICOMFileReader
::AnalyzeInputFiles()
{
this->SetGroup3DandT(true);
Superclass::AnalyzeInputFiles();
// collect output from superclass
size_t number_of_outputs = this->GetNumberOfOutputs();
if(number_of_outputs == 0)
{
MITK_ERROR << "Failed to parse input, retrieved 0 outputs from SeriesGDCMReader ";
}
MITK_INFO("diffusion.dicomreader") << "Retrieved " << number_of_outputs << "outputs.";
for( unsigned int outputidx = 0; outputidx < this->GetNumberOfOutputs(); outputidx++ )
{
DICOMImageBlockDescriptor block_0 = this->GetOutput(outputidx);
// collect vendor ID from the first output, first image
StringList inputFilename;
DICOMImageFrameInfo::Pointer frame_0 = block_0.GetImageFrameList().at(0);
inputFilename.push_back( frame_0->Filename );
mitk::DiffusionHeaderDICOMFileReader::Pointer headerReader;
bool isMosaic = false;
gdcm::Scanner gdcmScanner;
gdcm::Tag t_vendor(0x008, 0x0070);
gdcm::Tag t_imagetype(0x008, 0x008);
// add DICOM Tag for vendor
gdcmScanner.AddTag( t_vendor );
// add DICOM Tag for image type
gdcmScanner.AddTag( t_imagetype );
if( gdcmScanner.Scan( inputFilename ) )
{
// retrieve both vendor and image type
std::string vendor = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_vendor );
std::string image_type = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_imagetype );
MITK_INFO("diffusion.dicomreader") << "Output " << outputidx+1 << " Got vendor: " << vendor << " image type " << image_type;
// parse vendor tag
if( vendor.find("SIEMENS") != std::string::npos )
//&& image_type.find("DIFFUSION") != std::string::npos )
{
if( image_type.find("MOSAIC") != std::string::npos )
{
headerReader = mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::New();
isMosaic = true;
}
else
{
headerReader = mitk::DiffusionHeaderSiemensDICOMFileReader::New();
}
}
else if( vendor.find("GE") != std::string::npos )
{
headerReader = mitk::DiffusionHeaderGEDICOMFileReader::New();
}
else if( vendor.find("Philips") != std::string::npos )
{
headerReader = mitk::DiffusionHeaderPhilipsDICOMFileReader::New();
}
else
{
// unknown vendor
}
if( headerReader.IsNull() )
{
MITK_ERROR << "No header reader for given vendor. ";
continue;
}
}
else
{
continue;
}
bool canread = false;
// iterate over the threeD+t block
int numberOfTimesteps = block_0.GetIntProperty("timesteps", 1);
int framesPerTimestep = block_0.GetImageFrameList().size() / numberOfTimesteps;
for( int idx = 0; idx < numberOfTimesteps; idx++ )
{
int access_idx = idx * framesPerTimestep;
DICOMImageFrameInfo::Pointer frame = this->GetOutput( outputidx ).GetImageFrameList().at( access_idx );
canread = headerReader->ReadDiffusionHeader( frame->Filename );
}
if( canread )
{
// collect the information
mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType retrievedHeader = headerReader->GetHeaderInformation();
m_IsMosaicData.push_back(isMosaic);
m_OutputHeaderContainer.push_back( retrievedHeader );
m_OutputReaderContainer.push_back( headerReader );
}
}
this->SetNumberOfOutputs( this->m_OutputHeaderContainer.size() );
for( unsigned int outputidx = 0; outputidx < this->GetNumberOfOutputs(); outputidx++ )
{
// TODO : Analyze outputs + header information, i.e. for the loading confidence
MITK_INFO("diffusion.dicomreader") << "---- DICOM Analysis Report ---- :: Output " << outputidx+1 << " of " << this->GetNumberOfOutputs();
try{
PerformHeaderAnalysis( this->m_OutputHeaderContainer.at( outputidx) );
}
catch( const std::exception& se)
{
MITK_ERROR << "STD Exception " << se.what();
}
MITK_INFO("diffusion.dicomreader") << "===========================================";
}
}
bool mitk::DiffusionDICOMFileReader
::CanHandleFile(const std::string & /* filename */)
{
//FIXME :
return true;
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h
index 6cd30e36d7..e4da4b1849 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h
@@ -1,283 +1,281 @@
#ifndef MITKDIFFUSIONDICOMFILEREADERHELPER_H
#define MITKDIFFUSIONDICOMFILEREADERHELPER_H
-#include "mitkDiffusionImage.h"
-
#include "itkImageSeriesReader.h"
#include "itkVectorImage.h"
#include "itkImageRegionIteratorWithIndex.h"
namespace mitk
{
/**
* @brief The MosaicDescriptor struct is a help struct holding the necessary information for
* loading a mosaic DICOM file into an MITK file with correct geometry information
*/
struct MosaicDescriptor
{
unsigned int nimages;
bool slicenormalup;
itk::ImageBase<3>::SpacingType spacing;
itk::ImageBase<3>::DirectionType direction;
float origin[3];
};
class DiffusionDICOMFileReaderHelper
{
public:
typedef std::vector< std::string > StringContainer;
typedef std::vector< StringContainer > VolumeFileNamesContainer;
std::string GetOutputName(const VolumeFileNamesContainer& filenames)
{
typedef itk::Image< short, 3> InputImageType;
typedef itk::ImageSeriesReader< InputImageType > SeriesReaderType;
SeriesReaderType::Pointer probe_reader = SeriesReaderType::New();
probe_reader->SetFileNames( filenames.at(0) );
probe_reader->GenerateOutputInformation();
probe_reader->Update();
std::string seriesDescTag, seriesNumberTag, patientName;
SeriesReaderType::DictionaryArrayRawPointer inputDict = probe_reader->GetMetaDataDictionaryArray();
if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0008|103e", seriesDescTag ) )
seriesDescTag = "UNSeries";
if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0020|0011", seriesNumberTag ) )
seriesNumberTag = "00000";
if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0010|0010", patientName ) )
patientName = "UnknownName";
std::stringstream ss;
ss << seriesDescTag << "_" << seriesNumberTag << "_" << patientName;
return ss.str();
}
template< typename PixelType, unsigned int VecImageDimension>
typename itk::VectorImage< PixelType, VecImageDimension >::Pointer LoadToVector(
const VolumeFileNamesContainer& filenames
//const itk::ImageBase<3>::RegionType requestedRegion
)
{
typedef itk::Image< PixelType, 3> InputImageType;
typedef itk::ImageSeriesReader< InputImageType > SeriesReaderType;
typename SeriesReaderType::Pointer probe_reader = SeriesReaderType::New();
probe_reader->SetFileNames( filenames.at(0) );
probe_reader->GenerateOutputInformation();
const itk::ImageBase<3>::RegionType requestedRegion = probe_reader->GetOutput()->GetLargestPossibleRegion();
MITK_INFO << " --- Probe reader : \n" <<
" Retrieved LPR " << requestedRegion;
typedef itk::VectorImage< PixelType, 3 > VectorImageType;
typename VectorImageType::Pointer output_image = VectorImageType::New();
output_image->SetNumberOfComponentsPerPixel( filenames.size() );
output_image->SetSpacing( probe_reader->GetOutput()->GetSpacing() );
output_image->SetOrigin( probe_reader->GetOutput()->GetOrigin() );
output_image->SetDirection( probe_reader->GetOutput()->GetDirection() );
output_image->SetLargestPossibleRegion( probe_reader->GetOutput()->GetLargestPossibleRegion() );
output_image->SetBufferedRegion( requestedRegion );
output_image->Allocate();
itk::ImageRegionIterator< VectorImageType > vecIter(
output_image, requestedRegion );
VolumeFileNamesContainer::const_iterator volumesFileNamesIter = filenames.begin();
// iterate over the given volumes
unsigned int component = 0;
while( volumesFileNamesIter != filenames.end() )
{
MITK_INFO << " ======== Loading volume " << component+1 << " of " << filenames.size();
typename SeriesReaderType::Pointer volume_reader = SeriesReaderType::New();
volume_reader->SetFileNames( *volumesFileNamesIter );
try
{
volume_reader->UpdateLargestPossibleRegion();
}
catch( const itk::ExceptionObject &e)
{
mitkThrow() << " ITK Series reader failed : "<< e.what();
}
itk::ImageRegionConstIterator< InputImageType > iRCIter (
volume_reader->GetOutput(),
volume_reader->GetOutput()->GetLargestPossibleRegion() );
// transfer to vector image
iRCIter.GoToBegin();
vecIter.GoToBegin();
while( !iRCIter.IsAtEnd() )
{
typename VectorImageType::PixelType vector_pixel = vecIter.Get();
vector_pixel.SetElement( component, iRCIter.Get() );
vecIter.Set( vector_pixel );
++vecIter;
++iRCIter;
}
++volumesFileNamesIter;
component++;
}
return output_image;
}
/**
* Create the vector image for the resulting diffusion image from a mosaic DICOM image set,
* The method needs to be provided with the MosaicDescriptor struct to be able to compute the
* correct index and to set the geometry information of the image itself.
*/
template< typename PixelType, unsigned int VecImageDimension>
typename itk::VectorImage< PixelType, VecImageDimension >::Pointer LoadMosaicToVector(
const VolumeFileNamesContainer& filenames,
const MosaicDescriptor& mosaicInfo )
{
typedef itk::Image< PixelType, 3> MosaicImageType;
typedef itk::ImageFileReader< MosaicImageType > SingleImageReaderType;
// generate output
typedef itk::VectorImage< PixelType, 3 > VectorImageType;
VolumeFileNamesContainer::const_iterator volumesFileNamesIter = filenames.begin();
// probe the first file to retrieve the size of the 2d image
// we need this information to compute the index relation between mosaic and resulting 3d position
// but we need it only once
typename SingleImageReaderType::Pointer mosaic_probe = SingleImageReaderType::New();
mosaic_probe->SetFileName( (*volumesFileNamesIter).at(0) );
try
{
mosaic_probe->UpdateLargestPossibleRegion();
}
catch( const itk::ExceptionObject &e)
{
mitkThrow() << " ITK Image file reader failed : "<< e.what();
}
typename MosaicImageType::RegionType mosaic_lpr = mosaic_probe->GetOutput()->GetLargestPossibleRegion();
MITK_INFO << " == MOSAIC: " << mosaic_lpr;
itk::ImageBase<3>::SizeValueType images_per_row = ceil( sqrt( (float) mosaicInfo.nimages ) );
itk::ImageBase<3>::RegionType requestedRegion;
requestedRegion.SetSize( 0, mosaic_lpr.GetSize()[0]/images_per_row);
requestedRegion.SetSize( 1, mosaic_lpr.GetSize()[1]/images_per_row);
requestedRegion.SetSize( 2, mosaicInfo.nimages);
typename VectorImageType::Pointer output_image = VectorImageType::New();
output_image->SetNumberOfComponentsPerPixel( filenames.size() );
typename VectorImageType::DirectionType dmatrix;
dmatrix.SetIdentity();
std::vector<double> dirx = mosaic_probe->GetImageIO()->GetDirection(0);
std::vector<double> diry = mosaic_probe->GetImageIO()->GetDirection(1);
std::vector<double> dirz = mosaic_probe->GetImageIO()->GetDirection(2);
dmatrix.GetVnlMatrix().set_column( 0, &dirx[0] );
dmatrix.GetVnlMatrix().set_column( 1, &diry[0] );
dmatrix.GetVnlMatrix().set_column( 2, &dirz[0] );
/*
FIXME!!! The struct currently does not provide the geometry information
the loading works as required*/
output_image->SetSpacing( mosaicInfo.spacing );
output_image->SetOrigin( mosaic_probe->GetOutput()->GetOrigin() );
output_image->SetDirection( dmatrix );
output_image->SetLargestPossibleRegion( requestedRegion );
output_image->SetBufferedRegion( requestedRegion );
output_image->Allocate();
itk::ImageRegionIteratorWithIndex< VectorImageType > vecIter(
output_image, requestedRegion );
// hold the image sizes in an extra variable ( used very often )
typename MosaicImageType::SizeValueType dx = requestedRegion.GetSize()[0];
typename MosaicImageType::SizeValueType dy = requestedRegion.GetSize()[1];
// iterate over the given volumes
unsigned int component = 0;
while( volumesFileNamesIter != filenames.end() )
{
MITK_INFO << " ======== Loading volume " << component+1 << " of " << filenames.size();
typename SingleImageReaderType::Pointer mosaic_reader = SingleImageReaderType::New();
mosaic_reader->SetFileName( (*volumesFileNamesIter).at(0) );
try
{
mosaic_reader->UpdateLargestPossibleRegion();
}
catch( const itk::ExceptionObject &e)
{
mitkThrow() << " ITK Image file reader failed : "<< e.what();
}
typename MosaicImageType::Pointer current_mosaic = mosaic_reader->GetOutput();
vecIter.GoToBegin();
while( !vecIter.IsAtEnd() )
{
typename VectorImageType::PixelType vector_pixel = vecIter.Get();
typename VectorImageType::IndexType threeD_index = vecIter.GetIndex();
typename MosaicImageType::IndexType mosaic_index;
mosaic_index[2] = 0;
// first find the corresponding tile in the mosaic
// this is defined by the z-position of the vector (3D) image iterator
// in x : z_index % #images_in_grid
// in y : z_index / #images_in_grid
//
// the remaining is just computing the correct position in the mosaic, done by
// --------- index of (0,0,z) ----- + --- current 2d position ---
mosaic_index[0] = (threeD_index[2] % images_per_row) * dx + threeD_index[0];
mosaic_index[1] = (threeD_index[2] / images_per_row) * dy + threeD_index[1];
typename MosaicImageType::PixelType mosaic_pixel = current_mosaic->GetPixel( mosaic_index );
vector_pixel.SetElement( component, mosaic_pixel );
vecIter.Set( vector_pixel );
++vecIter;
}
++volumesFileNamesIter;
component++;
}
return output_image;
}
};
}
#endif // MITKDIFFUSIONDICOMFILEREADERHELPER_H
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h
index ed62cf3f1e..7e1f348f36 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h
@@ -1,117 +1,120 @@
/*===================================================================
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 MITKDIFFUSIONHEADERFILEREADER_H
#define MITKDIFFUSIONHEADERFILEREADER_H
#include <MitkDiffusionCoreExports.h>
+#include <mitkCommon.h>
#include <itkLightObject.h>
-#include "mitkDiffusionImage.h"
+#include <itkObjectFactory.h>
#include "gdcmScanner.h"
#include "gdcmReader.h"
+#include <vnl/vnl_vector_fixed.h>
+
namespace mitk
{
/**
* @brief The DiffusionImageHeaderInformation struct
*/
struct DiffusionImageDICOMHeaderInformation
{
DiffusionImageDICOMHeaderInformation()
: b_value(0),
baseline(false),
isotropic(false)
{
g_vector.fill(0);
}
void Print()
{
MITK_INFO << " DiffusionImageHeaderInformation : \n"
<< " : b value : " << b_value << "\n"
<< " : gradient : " << g_vector << "\n"
<< " : isotropic : " << isotropic << "\n --- \n";
}
unsigned int b_value;
vnl_vector_fixed< double, 3> g_vector;
bool baseline;
bool isotropic;
};
struct DiffusionImageMosaicDICOMHeaderInformation
: public DiffusionImageDICOMHeaderInformation
{
unsigned long n_images;
bool slicenormalup;
};
/**
* @class DiffusionHeaderDICOMFileReader
*
* @brief Abstract class for all vendor specific diffusion file header reader
*
* To provide a diffusion header reader for a new vendor, reimplement the \sa ReadDiffusionHeader method.
*/
class MitkDiffusionCore_EXPORT DiffusionHeaderDICOMFileReader
: public itk::LightObject
{
public:
typedef std::vector< DiffusionImageDICOMHeaderInformation > DICOMHeaderListType;
mitkClassMacro( DiffusionHeaderDICOMFileReader, itk::LightObject )
itkSimpleNewMacro( Self )
/**
* @brief IsDiffusionHeader Parse the given dicom file and collect the special diffusion image information
* @return
*/
virtual bool ReadDiffusionHeader( std::string ){ return false; }
DICOMHeaderListType GetHeaderInformation();
protected:
DiffusionHeaderDICOMFileReader();
virtual ~DiffusionHeaderDICOMFileReader();
DICOMHeaderListType m_HeaderInformationList;
};
/**
* @brief Retrieve the value of a gdcm tag to the given string
*
* @param tag the gdcm::Tag to be search for
* @param dataset a gdcm::DataSet to look into
* @param target a string to store the value of the given tag if found
* @param verbose make some output
*
* @return true if a string was found, false otherwise
*/
bool RevealBinaryTag(const gdcm::Tag tag, const gdcm::DataSet& dataset, std::string& target);
bool RevealBinaryTagC(const gdcm::Tag tag, const gdcm::DataSet& dataset, char* target_array );
} // end namespace mitk
#endif // MITKDIFFUSIONHEADERFILEREADER_H
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp
index 04f552619e..e90e9b77a4 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp
@@ -1,101 +1,102 @@
#include "mitkDiffusionHeaderGEDICOMFileReader.h"
+#include <vnl/vnl_math.h>
mitk::DiffusionHeaderGEDICOMFileReader
::DiffusionHeaderGEDICOMFileReader()
{
}
mitk::DiffusionHeaderGEDICOMFileReader
::~DiffusionHeaderGEDICOMFileReader()
{
}
bool mitk::DiffusionHeaderGEDICOMFileReader
::ReadDiffusionHeader(std::string filename)
{
gdcm::Reader gdcmReader;
gdcmReader.SetFileName( filename.c_str() );
gdcmReader.Read();
gdcm::Tag ge_bvalue_tag( 0x0043, 0x1039 );
gdcm::Tag ge_gradient_x( 0x0019, 0x10bb );
gdcm::Tag ge_gradient_y( 0x0019, 0x10bc );
gdcm::Tag ge_gradient_z( 0x0019, 0x10bd );
bool success = true;
DiffusionImageDICOMHeaderInformation header_info;
std::string ge_tagvalue_string;
char* pEnd;
// start with b-value
success = RevealBinaryTag( ge_bvalue_tag, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string );
// b value stored in the first bytes
// typical example: "1000\8\0\0" for bvalue=1000
// "40\8\0\0" for bvalue=40
// so we need to cut off the last 6 elements
const char* bval_string = ge_tagvalue_string.substr(0,ge_tagvalue_string.length()-6).c_str();
header_info.b_value = static_cast<unsigned int>(strtod( bval_string, &pEnd ));
// now retrieve the gradient direction
if(success &&
RevealBinaryTag( ge_gradient_x, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) )
{
header_info.g_vector[0] = strtod( ge_tagvalue_string.c_str(), &pEnd );
}
else
{
success = false;
}
if( success &&
RevealBinaryTag( ge_gradient_y, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) )
{
header_info.g_vector[1] = strtod( ge_tagvalue_string.c_str(), &pEnd );
}
else
{
success = false;
}
if( success &&
RevealBinaryTag( ge_gradient_z, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) )
{
header_info.g_vector[2] = strtod( ge_tagvalue_string.c_str(), &pEnd );
}
else
{
success = false;
}
if( success )
{
// Fix for (0,0,0) direction in IVIM datasets
if( header_info.b_value > 0 &&
header_info.g_vector.two_norm() < vnl_math::eps )
{
header_info.g_vector.fill(1);
header_info.g_vector.normalize();
header_info.isotropic = true;
}
// mark baseline
if( header_info.b_value == 0 )
header_info.baseline = true;
this->m_HeaderInformationList.push_back( header_info );
header_info.Print();
}
return success;
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp
index b037bbb3da..a9c8710521 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp
@@ -1,95 +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.
===================================================================*/
#include "mitkDiffusionHeaderPhilipsDICOMFileReader.h"
#include <sstream>
+#include <vnl/vnl_math.h>
+
mitk::DiffusionHeaderPhilipsDICOMFileReader::DiffusionHeaderPhilipsDICOMFileReader()
{
}
mitk::DiffusionHeaderPhilipsDICOMFileReader::~DiffusionHeaderPhilipsDICOMFileReader()
{
}
bool mitk::DiffusionHeaderPhilipsDICOMFileReader::ReadDiffusionHeader(std::string filename)
{
gdcm::Reader gdcmReader;
gdcmReader.SetFileName( filename.c_str() );
gdcmReader.Read();
gdcm::Tag philips_bvalue_tag( 0x2001, 0x1003 );
//gdcm::Tag philips_gradient_direction( 0x2001, 0x1004 );
DiffusionImageDICOMHeaderInformation header_info;
//std::string tagvalue_string;
//char* pEnd;
// reveal b-value
float bvalue = 0;
if( RevealBinaryTagC( philips_bvalue_tag, gdcmReader.GetFile().GetDataSet(), (char*) &bvalue) )
{
header_info.b_value = std::ceil( bvalue );
if( header_info.b_value == 0)
header_info.baseline = true;
}
else
{
MITK_WARN("diffusion.dicomreader.philips") << "No b-value found. Most probably no diffusion-weighted image.";
return false;
}
gdcm::Tag philips_new_bvalue_tag( 0x0018,0x9087 );
double dbvalue = 0;
if( RevealBinaryTagC( philips_new_bvalue_tag, gdcmReader.GetFile().GetDataSet(), (char*) &dbvalue) )
{
MITK_INFO("philips.dicom.diffusion.bvalue") << dbvalue;
}
if( header_info.baseline )
{
// no direction in unweighted images
header_info.g_vector.fill(0);
}
else
{
MITK_INFO("philips.dicom.diffusion.gradientdir") << "Parsing gradient direction.";
gdcm::Tag philips_gradient_direction_new( 0x0018, 0x9089 );
double gr_dir_arr[3] = {1,0,-1};
if( RevealBinaryTagC( philips_gradient_direction_new, gdcmReader.GetFile().GetDataSet(), (char*) &gr_dir_arr ) )
{
MITK_INFO("philips.dicom.diffusion.gradient") << "(" << gr_dir_arr[0] <<"," << gr_dir_arr[1] <<"," <<gr_dir_arr[2] <<")";
}
header_info.g_vector.copy_in( &gr_dir_arr[0] );
if( header_info.g_vector.two_norm() < vnl_math::eps )
{
header_info.g_vector.fill(1);
header_info.isotropic = true;
}
}
this->m_HeaderInformationList.push_back( header_info );
return true;
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp
index 043f5e94f2..d5e17f2e14 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp
@@ -1,42 +1,42 @@
#include "mitkDiffusionHeaderSiemensDICOMFileHelper.h"
#include <cstdlib>
#include <cmath>
mitk::SiemensDiffusionHeaderType mitk::GetHeaderType( std::string header )
{
// The CSA2 format begins with the string ‘SV10’, the CSA1 format does not.
if( header.find("SV10") != std::string::npos )
{
return mitk::SIEMENS_CSA2;
}
else
{
return mitk::SIEMENS_CSA1;
}
}
bool mitk::ParseInputString( std::string input, std::vector<double>& values, Siemens_Header_Format format_specs )
{
-
// TODO : Compute offset based on the format_specs, where does the 84 come from???
int offset = 84;
int vm = *(input.c_str() + format_specs.NameLength );
for (int k = 0; k < vm; k++)
{
int itemLength = *(input.c_str() + offset + 4);
int strideSize = static_cast<int> (ceil(static_cast<double>(itemLength)/4) * 4);
std::string valueString = input.substr( offset+16, itemLength );
double value = atof( valueString.c_str() );
values.push_back( value );
offset += 16+strideSize;
}
- return true;
+ // If there are no values it is invalid
+ return (values.size() > 0 );
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp
index 62416888ca..5cf26e5b8e 100644
--- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp
@@ -1,151 +1,152 @@
/*===================================================================
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 "mitkDiffusionHeaderSiemensDICOMFileReader.h"
#include "mitkDiffusionHeaderSiemensDICOMFileHelper.h"
#include "gdcmScanner.h"
#include "gdcmReader.h"
+#include <vnl/vnl_math.h>
/**
* @brief Extract b value from the siemens diffusion tag
*/
bool mitk::DiffusionHeaderSiemensDICOMFileReader
::ExtractSiemensDiffusionTagInformation( std::string tag_value, mitk::DiffusionImageDICOMHeaderInformation& values)
{
SiemensDiffusionHeaderType hformat = mitk::GetHeaderType( tag_value );
Siemens_Header_Format specs = this->m_SiemensFormatsCollection.at( hformat );
MITK_DEBUG << " Header format: " << hformat;
MITK_DEBUG << " :: Retrieving b value. ";
std::string::size_type tag_position =
tag_value.find( "B_value", 0 );
if( tag_position == std::string::npos )
{
MITK_ERROR << "No b value information found. ";
return false;
}
std::string value_string = tag_value.substr( tag_position, tag_value.size() - tag_position + 1 );
std::vector<double> value_array;
if( ParseInputString(value_string, value_array, specs) )
{
values.b_value = value_array.at(0);
}
else
{
MITK_INFO("diffusion.dicomreader.siemens") << "No b-value tag found. ";
return false;
}
// search for GradientDirectionInformation if the bvalue is not null
if( values.b_value > 0 )
{
std::string::size_type tag_position = tag_value.find( "DiffusionGradientDirection", 0 );
// Possibly it is a IVIM dataset, i.e. the gradient direction is not relevant
// and possibly either not set or set to zero
if( tag_position == std::string::npos )
{
MITK_WARN << "No gradient direction information, but non-zero b-value. Possibly an IVIM dataset. " << "\n"
<< "Setting gradient to (1,1,1).";
values.isotropic = true;
values.g_vector.fill(1);
return false;
}
value_array.clear();
std::string gradient_direction_str = tag_value.substr( tag_position, tag_value.size() - tag_position + 1 );
if( ParseInputString(gradient_direction_str, value_array, specs) )
{
if( value_array.size() != 3 )
{
MITK_ERROR << " Retrieved gradient information of length " << value_array.size();
return false;
}
for( unsigned int i=0; i<value_array.size(); i++)
{
values.g_vector[i] = value_array.at(i);
}
// Test for isotropic data (i.e. IVIM)
if( values.g_vector.two_norm() < vnl_math::eps )
{
values.g_vector.fill(1);
values.isotropic = true;
}
}
}
else
{
values.baseline = true;
}
return true;
}
mitk::DiffusionHeaderSiemensDICOMFileReader
::DiffusionHeaderSiemensDICOMFileReader()
{
Siemens_Header_Format Siemens_CSA1_Format( 64, sizeof(int32_t), 4, sizeof(int32_t), sizeof(int32_t) );
Siemens_Header_Format Siemens_CSA2_Format( 64, sizeof(int32_t), 4, sizeof(int32_t), sizeof(int32_t) );
m_SiemensFormatsCollection.push_back( Siemens_CSA1_Format );
m_SiemensFormatsCollection.push_back( Siemens_CSA2_Format );
}
mitk::DiffusionHeaderSiemensDICOMFileReader
::~DiffusionHeaderSiemensDICOMFileReader()
{
}
bool mitk::DiffusionHeaderSiemensDICOMFileReader
::ReadDiffusionHeader(std::string filename)
{
gdcm::Reader gdcmReader;
gdcmReader.SetFileName( filename.c_str() );
gdcmReader.Read();
MITK_INFO << " -- Analyzing: " << filename;
const gdcm::DataSet& dataset = gdcmReader.GetFile().GetDataSet();
const gdcm::Tag t_sie_diffusion( 0x0029,0x1010 );
//const gdcm::Tag t_sie_diffusion_vec( 0x0029,0x100e );
//const gdcm::Tag t_sie_diffusion2( 0x0029,0x100c );
std::string siemens_diffusionheader_str;
if( RevealBinaryTag( t_sie_diffusion, dataset, siemens_diffusionheader_str ) )
{
DiffusionImageDICOMHeaderInformation values;
this->ExtractSiemensDiffusionTagInformation( siemens_diffusionheader_str, values );
m_HeaderInformationList.push_back( values );
}
return true;
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp
index 2f909c5b73..d93af7cf2d 100644
--- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp
@@ -1,117 +1,118 @@
/*===================================================================
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 <mitkDiffusionPropertyHelper.h>
#include <vnl/algo/vnl_symmetric_eigensystem.h>
#include <vnl/vnl_inverse.h>
-template< typename TPixelType >
-mitk::DiffusionImageCorrectionFilter<TPixelType>::DiffusionImageCorrectionFilter()
+
+mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter()
{
}
-template< typename TPixelType >
-typename mitk::DiffusionImageCorrectionFilter<TPixelType>::TransformMatrixType
-mitk::DiffusionImageCorrectionFilter<TPixelType>
+
+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<TPixelType>
+
+void mitk::DiffusionImageCorrectionFilter
::CorrectDirections( const TransformsVectorType& transformations)
{
if( m_SourceImage.IsNull() )
{
mitkThrow() << " No diffusion image given! ";
}
- GradientDirectionContainerPointerType directions = m_SourceImage->GetDirections();
+ GradientDirectionContainerPointerType directions = static_cast<mitk::GradientDirectionsProperty*>( m_SourceImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
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 );
+ m_SourceImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( corrected_directions ) );
}
-template< typename TPixelType >
-void mitk::DiffusionImageCorrectionFilter<TPixelType>
+
+void mitk::DiffusionImageCorrectionFilter
::CorrectDirections( const TransformMatrixType& transformation)
{
if( m_SourceImage.IsNull() )
{
mitkThrow() << " No diffusion image given! ";
}
TransformsVectorType transfVec;
- for (unsigned int i=0; i< m_SourceImage->GetDirections()->Size();i++)
+ for (unsigned int i=0; i< static_cast<mitk::GradientDirectionsProperty*>( m_SourceImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->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 7d6c3b3e59..405cea700d 100644
--- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h
+++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h
@@ -1,97 +1,96 @@
/*===================================================================
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"
+#include "mitkImageSource.h"
+#include <MitkDiffusionCoreExports.h>
namespace mitk
{
/**
* @class DiffusionImageCorrectionFilter
*/
-template< typename TPixelType >
-class DiffusionImageCorrectionFilter
- : public DiffusionImageSource< TPixelType >
+class MitkDiffusionCore_EXPORT DiffusionImageCorrectionFilter
+ : public ImageSource
{
public:
/** class macros */
mitkClassMacro( DiffusionImageCorrectionFilter,
- DiffusionImageSource<TPixelType> )
+ ImageSource )
itkSimpleNewMacro(Self)
+ typedef short DiffusionPixelType;
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<TPixelType,3> ImageType;
+ typedef Superclass::OutputType DiffusionImageType;
+ typedef DiffusionImageType::Pointer DiffusionImageTypePointer;
+ typedef itk::VectorImage<DiffusionPixelType,3> 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
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake
index 4a2a9d22b0..9e87b8b605 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake
@@ -1,16 +1,15 @@
set(MODULE_TESTS
- mitkDiffusionImageEqualTest.cpp
mitkNonLocalMeansDenoisingTest.cpp
- mitkDiffusionImageEqualTest.cpp
+ mitkDiffusionPropertySerializerTest.cpp
)
set(MODULE_CUSTOM_TESTS
mitkPyramidImageRegistrationMethodTest.cpp
mitkDWHeadMotionCorrectionTest.cpp
mitkImageReconstructionTest.cpp
mitkConvertDWITypeTest.cpp
mitkExtractSingleShellTest.cpp
mitkNonLocalMeansDenoisingTest.cpp
mitkDiffusionDICOMFileReaderTest.cpp
)
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp
index c157efef62..101143f9d9 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp
@@ -1,109 +1,108 @@
/*===================================================================
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 "mitkTestingMacros.h"
#include "itkVectorImage.h"
#include "mitkBaseData.h"
#include "mitkBaseDataIOFactory.h"
#include "mitkImageWriter.h"
#include "itkImageFileWriter.h"
#include "itkB0ImageExtractionToSeparateImageFilter.h"
-#include "mitkDiffusionImage.h"
#include "mitkDiffusionCoreObjectFactory.h"
/** Documentation
* Test for factory registration
*/
int mitkB0ExtractionToSeparateImagesFilterTest(int argc , char* argv[])
{
// always start with this!
MITK_TEST_BEGIN("mitkB0ExtractionToSeparateImagesFilterTest");
MITK_TEST_CONDITION_REQUIRED(argc > 2, "Test image is specified. ");
typedef short DiffusionPixelType;
- typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType;
+ typedef mitk::Image DiffusionImageType;
typedef mitk::NrrdDiffusionImageReader< DiffusionPixelType> DiffusionNrrdReaderType;
RegisterDiffusionImagingObjectFactory();
std::string inputFileName( argv[1] );
std::vector<mitk::BaseData::Pointer> inputBaseDataVector = mitk::BaseDataIO::LoadBaseDataFromFile( inputFileName, "","",false);
MITK_TEST_CONDITION_REQUIRED( inputBaseDataVector.size() > 0, "BaseDataIO returned non-empty vector.");
mitk::BaseData::Pointer baseData = inputBaseDataVector.at(0);
MITK_TEST_CONDITION_REQUIRED( baseData.IsNotNull(), "BaseData is not null")
DiffusionImageType* vols = dynamic_cast< DiffusionImageType* >(baseData.GetPointer());
MITK_TEST_CONDITION_REQUIRED( vols != NULL, "Casting basedata to diffusion image successfull." );
// filter
typedef itk::B0ImageExtractionToSeparateImageFilter< short, short> FilterType;
typename FilterType::Pointer filter = FilterType::New();
MITK_TEST_CONDITION_REQUIRED(filter.IsNotNull(), "Filter instance created. ");
filter->SetInput(vols->GetVectorImage());
filter->SetDirections(vols->GetDirections());
filter->Update();
// output
mitk::Image::Pointer mitkImage = mitk::Image::New();
MITK_TEST_CONDITION_REQUIRED( mitkImage.IsNotNull(), "mitkImage not null." );
mitkImage->InitializeByItk( filter->GetOutput() );
MITK_TEST_CONDITION_REQUIRED( mitkImage->GetDimension()==4, "Output image is a 4D image.");
mitkImage->SetImportChannel( filter->GetOutput()->GetBufferPointer() );
typedef itk::ImageFileWriter< FilterType::OutputImageType > itkImageWriterType;
typename itkImageWriterType::Pointer itkWriter = itkImageWriterType::New();
itkWriter->SetFileName( argv[2] );
itkWriter->SetInput( filter->GetOutput() );
try
{
itkWriter->Update();
}
catch(itk::ExceptionObject &e)
{
MITK_ERROR << "Catched exception from image writer. " << e.what();
}
/*
// write output
mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New();
MITK_TEST_CONDITION_REQUIRED( writer.IsNotNull(), "Writer instance created. ");
writer->SetInput( mitkImage );
writer->SetExtension(".nrrd");
writer->SetFileName( "/localdata/hering/_Images/TestB0Extraction" );
writer->Update();
*/
// always end with this!
MITK_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp
index aeae0f8536..4e6449576a 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp
@@ -1,51 +1,46 @@
/*===================================================================
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 "mitkTestingMacros.h"
#include "mitkIOUtil.h"
#include "mitkDWIHeadMotionCorrectionFilter.h"
-typedef short DiffusionPixelType;
-typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType;
-
/**
* @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter
*
* @param argv : Input and Output image full path
*/
int mitkConvertDWITypeTest( int argc, char* argv[] )
{
MITK_TEST_BEGIN("mitkConvertDWITypeTest");
MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output.");
mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] );
- DiffusionImageType* dwimage =
- static_cast<DiffusionImageType*>( inputImage.GetPointer() );
try
{
- mitk::IOUtil::Save(dwimage, argv[2]);
+ mitk::IOUtil::Save(inputImage, argv[2]);
}
catch( const itk::ExceptionObject& e)
{
- MITK_ERROR << "Catched exception: " << e.what();
+ MITK_ERROR << "Caught exception: " << e.what();
mitkThrow() << "Failed with exception from subprocess!";
}
MITK_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp
index 7805c72b3c..e19fc61f9c 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp
@@ -1,60 +1,60 @@
/*===================================================================
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 "mitkTestingMacros.h"
#include "mitkIOUtil.h"
-#include "mitkDWIHeadMotionCorrectionFilter.h"
+#include <mitkDWIHeadMotionCorrectionFilter.h>
typedef short DiffusionPixelType;
-typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType;
+typedef mitk::Image DiffusionImageType;
/**
* @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter
*
* @param argv : Input and Output image full path
*/
int mitkDWHeadMotionCorrectionTest( int argc, char* argv[] )
{
MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest");
MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output.");
// itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1);
mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] );
DiffusionImageType* dwimage =
static_cast<DiffusionImageType*>( inputImage.GetPointer() );
- mitk::DWIHeadMotionCorrectionFilter<DiffusionPixelType>::Pointer corrfilter =
- mitk::DWIHeadMotionCorrectionFilter<DiffusionPixelType>::New();
+ mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter =
+ mitk::DWIHeadMotionCorrectionFilter::New();
corrfilter->SetInput( dwimage );
corrfilter->SetAverageUnweighted(false);
corrfilter->Update();
try
{
mitk::IOUtil::SaveBaseData(corrfilter->GetOutput(), argv[2]);
}
catch( const itk::ExceptionObject& e)
{
MITK_ERROR << "Catched exception: " << e.what();
mitkThrow() << "Failed with exception from subprocess!";
}
MITK_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp
index 1b18949742..bfdd73c944 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp
@@ -1,115 +1,113 @@
/*===================================================================
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 "mitkDiffusionDICOMFileReader.h"
#include "mitkDiffusionDICOMFileReaderTestHelper.h"
#include "mitkDICOMTagBasedSorter.h"
#include "mitkDICOMSortByTag.h"
#include <mitkIOUtil.h>
#include "mitkTestingMacros.h"
using mitk::DICOMTag;
int mitkDiffusionDICOMFileReaderTest(int argc, char* argv[])
{
MITK_TEST_BEGIN("mitkDiffusionDICOMFileReaderTest");
mitk::DiffusionDICOMFileReader::Pointer gdcmReader = mitk::DiffusionDICOMFileReader::New();
MITK_TEST_CONDITION_REQUIRED(gdcmReader.IsNotNull(), "DICOMITKSeriesGDCMReader can be instantiated.");
std::string output_filename = "/tmp/dicom_out.dwi";
if( argc > 3)
{
mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc-1,argv );
output_filename = std::string( argv[argc-1] );
}
else
{
mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc,argv );
}
// check the Set/GetInput function
mitk::DICOMFileReaderTestHelper::TestInputFilenames( gdcmReader );
MITK_INFO << "Test input filenanems";
// check that output is a good reproduction of input (no duplicates, no new elements)
mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader );
MITK_INFO << "Test output";
// repeat test with some more realistic sorting
gdcmReader = mitk::DiffusionDICOMFileReader::New(); // this also tests destruction
mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New();
// Use tags as in Qmitk
// all the things that split by tag in DicomSeriesReader
tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0010) ); // Number of Rows
tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0011) ); // Number of Columns
tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0030) ); // Pixel Spacing
tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing
tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code)
// TODO handle as real vectors! cluster with configurable errors!
//tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x000e) ); // Series Instance UID
//tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0010) );
tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x0050) ); // Slice Thickness
tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0008) ); // Number of Frames
tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0052) ); // Frame of Reference UID
// gdcmReader->AddSortingElement( tagSorter );
//mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader );
mitk::DICOMSortCriterion::ConstPointer sorting =
mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0013), // instance number
mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0012), // aqcuisition number
mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0032), // aqcuisition time
mitk::DICOMSortByTag::New( DICOMTag(0x0018, 0x1060), // trigger time
mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly)
).GetPointer()
).GetPointer()
).GetPointer()
).GetPointer()
).GetPointer();
tagSorter->SetSortCriterion( sorting );
MITK_INFO << "Created sort";
gdcmReader->AddSortingElement( tagSorter );
mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader );
MITK_INFO << "Created sort";
//gdcmReader->PrintOutputs(std::cout, true);
// really load images
//mitk::DICOMFileReaderTestHelper::TestMitkImagesAreLoaded( gdcmReader );
gdcmReader->LoadImages();
mitk::Image::Pointer loaded_image = gdcmReader->GetOutput(0).GetMitkImage();
- mitk::DiffusionImage<short>::Pointer d_img = static_cast<mitk::DiffusionImage<short>*>( loaded_image.GetPointer() );
-
try
{
- mitk::IOUtil::Save(d_img, output_filename.c_str());
+ mitk::IOUtil::Save(loaded_image, output_filename.c_str());
}
catch( const itk::ExceptionObject& e)
{
MITK_TEST_FAILED_MSG( << "Writer failed : " << e.what() );
}
MITK_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp
new file mode 100644
index 0000000000..436646692c
--- /dev/null
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp
@@ -0,0 +1,184 @@
+/*===================================================================
+
+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 "mitkIOUtil.h"
+#include "mitkTestingMacros.h"
+#include "mitkTestFixture.h"
+
+#include <mitkBasePropertySerializer.h>
+
+#include <mitkBValueMapProperty.h>
+
+#include <mitkGradientDirectionsProperty.h>
+#include <mitkMeasurementFrameProperty.h>
+
+#include <mitkDiffusionPropertyHelper.h>
+
+class mitkDiffusionPropertySerializerTestSuite : public mitk::TestFixture
+{
+
+ CPPUNIT_TEST_SUITE(mitkDiffusionPropertySerializerTestSuite);
+ MITK_TEST(Equal_SerializeandDeserialize_ReturnsTrue);
+
+ //MITK_TEST(Equal_DifferentChannels_ReturnFalse);
+
+
+ CPPUNIT_TEST_SUITE_END();
+
+private:
+
+ /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/
+ mitk::PropertyList::Pointer propList; //represet image propertylist
+ mitk::BValueMapProperty::Pointer bvaluemap_prop;
+ mitk::GradientDirectionsProperty::Pointer gradientdirection_prop;
+ mitk::MeasurementFrameProperty::Pointer measurementframe_prop;
+
+public:
+
+ /**
+* @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called).
+*/
+ void setUp()
+ {
+
+ propList = mitk::PropertyList::New();
+
+ mitk::BValueMapProperty::BValueMap map;
+ std::vector<unsigned int> indices1;
+ indices1.push_back(1);
+ indices1.push_back(2);
+ indices1.push_back(3);
+ indices1.push_back(4);
+
+ map[0] = indices1;
+ std::vector<unsigned int> indices2;
+ indices2.push_back(4);
+ indices2.push_back(3);
+ indices2.push_back(2);
+ indices2.push_back(1);
+
+ map[1000] = indices2;
+ bvaluemap_prop = mitk::BValueMapProperty::New(map).GetPointer();
+ propList->SetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str(),bvaluemap_prop);
+
+
+ mitk::GradientDirectionsProperty::GradientDirectionsContainerType::Pointer gdc;
+ gdc = mitk::GradientDirectionsProperty::GradientDirectionsContainerType::New();
+
+ double a[3] = {3.0,4.0,1.4};
+ vnl_vector_fixed<double,3> vec1;
+ vec1.set(a);
+ gdc->push_back(vec1);
+
+ double b[3] = {1.0,5.0,123.4};
+ vnl_vector_fixed<double,3> vec2;
+ vec2.set(b);
+ gdc->push_back(vec2);
+
+ double c[3] = {13.0,84.02,13.4};
+ vnl_vector_fixed<double,3> vec3;
+ vec3.set(c);
+ gdc->push_back(vec3);
+
+
+ gradientdirection_prop = mitk::GradientDirectionsProperty::New(gdc.GetPointer()).GetPointer();
+ propList->SetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), gradientdirection_prop);
+
+ mitk::MeasurementFrameProperty::MeasurementFrameType mft;
+
+ double row0[3] = {1,0,0};
+ double row1[3] = {0,1,0};
+ double row2[3] = {0,0,1};
+
+ mft.set_row(0,row0);
+ mft.set_row(1,row1);
+ mft.set_row(2,row2);
+
+ measurementframe_prop = mitk::MeasurementFrameProperty::New(mft).GetPointer();
+ propList->SetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), measurementframe_prop);
+ }
+
+ void tearDown()
+ {
+
+ }
+
+ void Equal_SerializeandDeserialize_ReturnsTrue()
+ {
+
+ assert(propList);
+
+ /* try to serialize each property in the list, then deserialize again and check for equality */
+ for (mitk::PropertyList::PropertyMap::const_iterator it = propList->GetMap()->begin(); it != propList->GetMap()->end(); ++it)
+ {
+ const mitk::BaseProperty* prop = it->second;
+ // construct name of serializer class
+ std::string serializername = std::string(prop->GetNameOfClass()) + "Serializer";
+ std::list<itk::LightObject::Pointer> allSerializers = itk::ObjectFactoryBase::CreateAllInstance(serializername.c_str());
+ MITK_TEST_CONDITION(allSerializers.size() > 0, std::string("Creating serializers for ") + serializername);
+ if (allSerializers.size() == 0)
+ {
+ MITK_TEST_OUTPUT( << "serialization not possible, skipping " << prop->GetNameOfClass());
+ continue;
+ }
+ if (allSerializers.size() > 1)
+ {
+ MITK_TEST_OUTPUT (<< "Warning: " << allSerializers.size() << " serializers found for " << prop->GetNameOfClass() << "testing only the first one.");
+ }
+ mitk::BasePropertySerializer* serializer = dynamic_cast<mitk::BasePropertySerializer*>( allSerializers.begin()->GetPointer());
+ MITK_TEST_CONDITION(serializer != NULL, serializername + std::string(" is valid"));
+ if (serializer != NULL)
+ {
+ serializer->SetProperty(prop);
+ TiXmlElement* valueelement = NULL;
+ try
+ {
+ valueelement = serializer->Serialize();
+// TiXmlPrinter p;
+// valueelement->Accept(&p);
+// MITK_INFO << p.CStr();
+ }
+ catch (...)
+ {
+ }
+ MITK_TEST_CONDITION(valueelement != NULL, std::string("Serialize property with ") + serializername);
+
+ if (valueelement == NULL)
+ {
+ MITK_TEST_OUTPUT( << "serialization failed, skipping deserialization");
+ continue;
+ }
+
+ mitk::BaseProperty::Pointer deserializedProp = serializer->Deserialize( valueelement );
+ MITK_TEST_CONDITION(deserializedProp.IsNotNull(), "serializer created valid property");
+ if (deserializedProp.IsNotNull())
+ {
+ MITK_TEST_CONDITION(*(deserializedProp.GetPointer()) == *prop, "deserialized property equals initial property for type " << prop->GetNameOfClass());
+ }
+
+ }
+ else
+ {
+ MITK_TEST_OUTPUT( << "created serializer object is of class " << allSerializers.begin()->GetPointer()->GetNameOfClass())
+ }
+ } // for all properties
+
+ }
+
+
+};
+
+MITK_TEST_SUITE_REGISTRATION(mitkDiffusionPropertySerializer)
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp
index f2b1ae2338..43686ea03e 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp
@@ -1,96 +1,103 @@
/*===================================================================
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 "mitkTestingMacros.h"
#include "mitkIOUtil.h"
#include <itkElectrostaticRepulsionDiffusionGradientReductionFilter.h>
#include "mitkDWIHeadMotionCorrectionFilter.h"
+#include <mitkDiffusionPropertyHelper.h>
+#include <mitkImageCast.h>
+#include <mitkProperties.h>
typedef short DiffusionPixelType;
-typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType;
int mitkExtractSingleShellTest( int argc, char* argv[] )
{
MITK_TEST_BEGIN("mitkExtractSingleShellTest");
MITK_TEST_CONDITION_REQUIRED( argc > 3, "Specify input and output and the shell to be extracted");
/*
1. Get input data
*/
- mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] );
- DiffusionImageType* dwimage =
- static_cast<DiffusionImageType*>( inputImage.GetPointer() );
+ mitk::Image::Pointer dwimage = mitk::IOUtil::LoadImage( argv[1] );
- MITK_TEST_CONDITION_REQUIRED( dwimage != NULL, "Input is a dw-image");
+ mitk::GradientDirectionsProperty::Pointer gradientsProperty = static_cast<mitk::GradientDirectionsProperty *>( dwimage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() );
+
+ MITK_TEST_CONDITION_REQUIRED( gradientsProperty.IsNotNull(), "Input is a dw-image");
unsigned int extract_value = 0;
std::istringstream input(argv[3]);
input >> extract_value;
typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter<DiffusionPixelType, DiffusionPixelType> FilterType;
- typedef DiffusionImageType::BValueMap BValueMap;
+ typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMap;
// GetShellSelection from GUI
BValueMap shellSelectionMap;
- BValueMap originalShellMap = dwimage->GetBValueMap();
+ BValueMap originalShellMap = static_cast<mitk::BValueMapProperty*>(dwimage->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap();
std::vector<unsigned int> newNumGradientDirections;
shellSelectionMap[extract_value] = originalShellMap[extract_value];
newNumGradientDirections.push_back( originalShellMap[extract_value].size() ) ;
- DiffusionImageType::GradientDirectionContainerType::Pointer gradientContainer = dwimage->GetDirections();
+ itk::VectorImage< short, 3 >::Pointer itkVectorImagePointer = itk::VectorImage< short, 3 >::New();
+ mitk::CastToItkImage(dwimage, itkVectorImagePointer);
+ itk::VectorImage< short, 3 > *vectorImage = itkVectorImagePointer.GetPointer();
+
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainer = static_cast<mitk::GradientDirectionsProperty*>( dwimage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer();
FilterType::Pointer filter = FilterType::New();
- filter->SetInput(dwimage->GetVectorImage());
+ filter->SetInput(vectorImage);
filter->SetOriginalGradientDirections(gradientContainer);
filter->SetNumGradientDirections(newNumGradientDirections);
filter->SetOriginalBValueMap(originalShellMap);
filter->SetShellSelectionBValueMap(shellSelectionMap);
try
{
filter->Update();
}
catch( const itk::ExceptionObject& e)
{
MITK_TEST_FAILED_MSG( << "Failed due to ITK exception: " << e.what() );
}
- DiffusionImageType::Pointer image = DiffusionImageType::New();
- image->SetVectorImage( filter->GetOutput() );
- image->SetReferenceBValue(dwimage->GetReferenceBValue());
- image->SetDirections(filter->GetGradientDirections());
- image->SetMeasurementFrame(dwimage->GetMeasurementFrame());
- image->InitializeFromVectorImage();
+ mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetGradientDirections() ) );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast<mitk::FloatProperty*>(dwimage->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) );
+ outImage->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast<mitk::MeasurementFrameProperty*>(dwimage->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) );
+ mitk::DiffusionPropertyHelper propertyHelper( outImage );
+ propertyHelper.InitializeImage();
+
/*
* 3. Write output data
**/
try
{
- mitk::IOUtil::Save(image, argv[2]);
+ mitk::IOUtil::Save(outImage, argv[2]);
}
catch( const itk::ExceptionObject& e)
{
MITK_ERROR << "Catched exception: " << e.what();
mitkThrow() << "Failed with exception from subprocess!";
}
MITK_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp
index e1f16dcca7..33bfd4cd53 100755
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp
@@ -1,152 +1,158 @@
/*===================================================================
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 <mitkImageCast.h>
#include <mitkTensorImage.h>
#include <mitkQBallImage.h>
-#include <mitkDiffusionImage.h>
#include <mitkIOUtil.h>
#include <mitkBaseDataIOFactory.h>
#include <itkDiffusionTensor3D.h>
#include <mitkTestingMacros.h>
#include <itkDiffusionTensor3DReconstructionImageFilter.h>
#include <itkDiffusionTensor3D.h>
#include <itkDiffusionQballReconstructionImageFilter.h>
#include <itkAnalyticalDiffusionQballReconstructionImageFilter.h>
+#include <mitkImage.h>
+#include <mitkDiffusionPropertyHelper.h>
using namespace std;
int mitkImageReconstructionTest(int argc, char* argv[])
{
MITK_TEST_BEGIN("mitkImageReconstructionTest");
MITK_TEST_CONDITION_REQUIRED(argc>1,"check for input data")
try
{
- mitk::DiffusionImage<short>::Pointer dwi = dynamic_cast<mitk::DiffusionImage<short>*>(mitk::IOUtil::LoadDataNode(argv[1])->GetData());
+ mitk::Image::Pointer dwi = dynamic_cast<mitk::Image*>(mitk::IOUtil::LoadDataNode(argv[1])->GetData());
+ itk::VectorImage<short,3>::Pointer itkVectorImagePointer = itk::VectorImage<short,3>::New();
+ mitk::CastToItkImage(dwi, itkVectorImagePointer);
+
+ float b_value = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi );
+ mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradients = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi);
{
MITK_INFO << "Tensor reconstruction " << argv[2];
mitk::TensorImage::Pointer tensorImage = dynamic_cast<mitk::TensorImage*>(mitk::IOUtil::LoadDataNode(argv[2])->GetData());
typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->Update();
mitk::TensorImage::Pointer testImage = mitk::TensorImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *tensorImage, 0.0001, true), "tensor reconstruction test.");
}
{
MITK_INFO << "Numerical Q-ball reconstruction " << argv[3];
mitk::QBallImage::Pointer qballImage = dynamic_cast<mitk::QBallImage*>(mitk::IOUtil::LoadDataNode(argv[3])->GetData());
typedef itk::DiffusionQballReconstructionImageFilter<short, short, float, QBALL_ODFSIZE> QballReconstructionImageFilterType;
QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD);
filter->Update();
mitk::QBallImage::Pointer testImage = mitk::QBallImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Numerical Q-ball reconstruction test.");
}
{
MITK_INFO << "Standard Q-ball reconstruction " << argv[4];
mitk::QBallImage::Pointer qballImage = dynamic_cast<mitk::QBallImage*>(mitk::IOUtil::LoadDataNode(argv[4])->GetData());
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->SetLambda(0.006);
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
filter->Update();
mitk::QBallImage::Pointer testImage = mitk::QBallImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Standard Q-ball reconstruction test.");
}
{
MITK_INFO << "CSA Q-ball reconstruction " << argv[5];
mitk::QBallImage::Pointer qballImage = dynamic_cast<mitk::QBallImage*>(mitk::IOUtil::LoadDataNode(argv[5])->GetData());
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->SetLambda(0.006);
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
filter->Update();
mitk::QBallImage::Pointer testImage = mitk::QBallImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "CSA Q-ball reconstruction test.");
}
{
MITK_INFO << "ADC profile reconstruction " << argv[6];
mitk::QBallImage::Pointer qballImage = dynamic_cast<mitk::QBallImage*>(mitk::IOUtil::LoadDataNode(argv[6])->GetData());
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->SetLambda(0.006);
filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY);
filter->Update();
mitk::QBallImage::Pointer testImage = mitk::QBallImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "ADC profile reconstruction test.");
}
{
MITK_INFO << "Raw signal modeling " << argv[7];
mitk::QBallImage::Pointer qballImage = dynamic_cast<mitk::QBallImage*>(mitk::IOUtil::LoadDataNode(argv[7])->GetData());
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter<short,short,float,4,QBALL_ODFSIZE> FilterType;
FilterType::Pointer filter = FilterType::New();
- filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() );
- filter->SetBValue(dwi->GetReferenceBValue());
+ filter->SetGradientImage( gradients, itkVectorImagePointer );
+ filter->SetBValue( b_value );
filter->SetLambda(0.006);
filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL);
filter->Update();
mitk::QBallImage::Pointer testImage = mitk::QBallImage::New();
testImage->InitializeByItk( filter->GetOutput() );
testImage->SetVolume( filter->GetOutput()->GetBufferPointer() );
MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Raw signal modeling test.");
}
}
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_TEST_END();
}
diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp
index 680b5c6295..62a7941cb4 100644
--- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp
+++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp
@@ -1,176 +1,169 @@
/*===================================================================
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 "mitkDiffusionImage.h"
#include "mitkIOUtil.h"
#include "mitkTestingMacros.h"
#include "mitkTestFixture.h"
#include "itkNonLocalMeansDenoisingFilter.h"
+#include "mitkGradientDirectionsProperty.h"
+#include "mitkITKImageImport.h"
+#include <mitkImageCast.h>
class mitkNonLocalMeansDenoisingTestSuite : public mitk::TestFixture
{
CPPUNIT_TEST_SUITE(mitkNonLocalMeansDenoisingTestSuite);
MITK_TEST(Denoise_NLMg_shouldReturnTrue);
MITK_TEST(Denoise_NLMr_shouldReturnTrue);
MITK_TEST(Denoise_NLMv_shouldReturnTrue);
MITK_TEST(Denoise_NLMvr_shouldReturnTrue);
CPPUNIT_TEST_SUITE_END();
private:
+ typedef itk::VectorImage<short,3> VectorImagetType;
+
/** Members used inside the different (sub-)tests. All members are initialized via setUp().*/
- mitk::DiffusionImage<short>::Pointer m_Image;
- mitk::DiffusionImage<short>::Pointer m_ReferenceImage;
- mitk::DiffusionImage<short>::Pointer m_DenoisedImage;
+ mitk::Image::Pointer m_Image;
+ mitk::Image::Pointer m_ReferenceImage;
+ mitk::Image::Pointer m_DenoisedImage;
itk::Image<short, 3>::Pointer m_ImageMask;
itk::NonLocalMeansDenoisingFilter<short>::Pointer m_DenoisingFilter;
public:
/**
* @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called).
*/
void setUp()
{
//generate test images
std::string imagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi.dwi");
-// std::string maskPath = GetTestDataFilePath("DiffusionImaging/Denoising/denoising_mask.nrrd");
- m_Image = static_cast<mitk::DiffusionImage<short>*>( mitk::IOUtil::LoadImage(imagePath).GetPointer());
-// mitk::Image::Pointer imageMask = static_cast<mitk::Image*>( mitk::IOUtil::LoadImage(maskPath).GetPointer());
-// mitk::CastToItkImage(imageMask, m_ImageMask);
+ m_Image = mitk::IOUtil::LoadImage(imagePath);
m_ReferenceImage = NULL;
- m_DenoisedImage = mitk::DiffusionImage<short>::New();
+ m_DenoisedImage = mitk::Image::New();
//initialise Filter
m_DenoisingFilter = itk::NonLocalMeansDenoisingFilter<short>::New();
- m_DenoisingFilter->SetInputImage(m_Image->GetVectorImage());
-// m_DenoisingFilter->SetInputMask(m_ImageMask);
+ VectorImagetType::Pointer vectorImage;
+ mitk::CastToItkImage(m_Image,vectorImage);
+ m_DenoisingFilter->SetInputImage(vectorImage);
m_DenoisingFilter->SetNumberOfThreads(1);
m_DenoisingFilter->SetComparisonRadius(1);
m_DenoisingFilter->SetSearchRadius(1);
m_DenoisingFilter->SetVariance(500);
}
void tearDown()
{
m_Image = NULL;
m_ImageMask = NULL;
m_ReferenceImage = NULL;
m_DenoisingFilter = NULL;
m_DenoisedImage = NULL;
}
void Denoise_NLMg_shouldReturnTrue()
{
std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMg.dwi");
- m_ReferenceImage = static_cast<mitk::DiffusionImage<short>*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer());
+ m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath);
m_DenoisingFilter->SetUseRicianAdaption(false);
m_DenoisingFilter->SetUseJointInformation(false);
try
{
m_DenoisingFilter->Update();
}
catch(std::exception& e)
{
MITK_ERROR << e.what();
}
- m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput());
- m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue());
- m_DenoisedImage->SetDirections(m_Image->GetDirections());
- m_DenoisedImage->InitializeFromVectorImage();
+ mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage);
+ m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone());
MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMg should always return the same result.");
}
void Denoise_NLMr_shouldReturnTrue()
{
std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMr.dwi");
- m_ReferenceImage = static_cast<mitk::DiffusionImage<short>*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer());
+ m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath);
m_DenoisingFilter->SetUseRicianAdaption(true);
m_DenoisingFilter->SetUseJointInformation(false);
try
{
m_DenoisingFilter->Update();
}
catch(std::exception& e)
{
MITK_ERROR << e.what();
}
- m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput());
- m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue());
- m_DenoisedImage->SetDirections(m_Image->GetDirections());
- m_DenoisedImage->InitializeFromVectorImage();
+ mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage);
+ m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone());
MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMr should always return the same result.");
}
void Denoise_NLMv_shouldReturnTrue()
{
std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMv.dwi");
- m_ReferenceImage = static_cast<mitk::DiffusionImage<short>*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer());
-
+ m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath);
m_DenoisingFilter->SetUseRicianAdaption(false);
m_DenoisingFilter->SetUseJointInformation(true);
try
{
m_DenoisingFilter->Update();
}
catch(std::exception& e)
{
MITK_ERROR << e.what();
}
- m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput());
- m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue());
- m_DenoisedImage->SetDirections(m_Image->GetDirections());
- m_DenoisedImage->InitializeFromVectorImage();
+ mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage);
+ m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone());
MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMv should always return the same result.");
}
void Denoise_NLMvr_shouldReturnTrue()
{
std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMvr.dwi");
- m_ReferenceImage = static_cast<mitk::DiffusionImage<short>*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer());
+ m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath);
m_DenoisingFilter->SetUseRicianAdaption(true);
m_DenoisingFilter->SetUseJointInformation(true);
try
{
m_DenoisingFilter->Update();
}
catch(std::exception& e)
{
MITK_ERROR << e.what();
}
- m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput());
- m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue());
- m_DenoisedImage->SetDirections(m_Image->GetDirections());
- m_DenoisedImage->InitializeFromVectorImage();
+ mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage);
+ m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone());
MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMvr should always return the same result.");
}
};
MITK_TEST_SUITE_REGISTRATION(mitkNonLocalMeansDenoising)
diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake
index d5ad2057eb..dae157948b 100644
--- a/Modules/DiffusionImaging/DiffusionCore/files.cmake
+++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake
@@ -1,150 +1,149 @@
set(CPP_FILES
# DicomImport
DicomImport/mitkDicomDiffusionImageReader.cpp
# DicomImport/mitkGroupDiffusionHeadersFilter.cpp
DicomImport/mitkDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp
DicomImport/mitkDiffusionDICOMFileReader.cpp
DicomImport/mitkDiffusionHeaderDICOMFileReader.cpp
DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp
DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp
DicomImport/mitkDiffusionHeaderSiemensMosaicDICOMFileReader.cpp
DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp
DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp
# DataStructures -> DWI
IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp
- IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp
-
- IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp
IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp
# Properties
IODataStructures/Properties/mitkBValueMapProperty.cpp
IODataStructures/Properties/mitkGradientDirectionsProperty.cpp
IODataStructures/Properties/mitkMeasurementFrameProperty.cpp
IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp
IODataStructures/Properties/mitkNodePredicateIsDWI.cpp
# Serializer
IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp
IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp
IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp
# DataStructures -> QBall
IODataStructures/QBallImages/mitkQBallImageSource.cpp
IODataStructures/QBallImages/mitkQBallImage.cpp
# DataStructures -> Tensor
IODataStructures/TensorImages/mitkTensorImage.cpp
Rendering/vtkMaskedProgrammableGlyphFilter.cpp
Rendering/mitkVectorImageVtkGlyphMapper3D.cpp
Rendering/vtkOdfSource.cxx
Rendering/vtkThickPlane.cxx
Rendering/mitkOdfNormalizationMethodProperty.cpp
Rendering/mitkOdfScaleByProperty.cpp
# Algorithms
Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp
Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp
Algorithms/itkDwiGradientLengthCorrectionFilter.cpp
Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h
# Registration Algorithms & Co.
Algorithms/Registration/mitkRegistrationWrapper.cpp
Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp
# Algorithms/Registration/mitkRegistrationMethodITK4.cpp
+ Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp
# MultishellProcessing
Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.cpp
Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.cpp
# Function Collection
mitkDiffusionFunctionCollection.cpp
)
set(H_FILES
# function Collection
mitkDiffusionFunctionCollection.h
# Rendering
- Rendering/mitkDiffusionImageMapper.h
Rendering/mitkOdfVtkMapper2D.h
# Reconstruction
Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h
Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h
Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h
Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h
Algorithms/Reconstruction/itkPointShell.h
Algorithms/Reconstruction/itkOrientationDistributionFunction.h
Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h
# MultishellProcessing
Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h
Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h
Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h
- # IO Datastructures
- IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h
+ # Properties
+ IODataStructures/Properties/mitkBValueMapProperty.h
+ IODataStructures/Properties/mitkGradientDirectionsProperty.h
+ IODataStructures/Properties/mitkMeasurementFrameProperty.h
+ IODataStructures/Properties/mitkDiffusionPropertyHelper.h
# Algorithms
Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h
Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h
Algorithms/itkTensorDerivedMeasurementsFilter.h
Algorithms/itkBrainMaskExtractionImageFilter.h
Algorithms/itkB0ImageExtractionImageFilter.h
Algorithms/itkB0ImageExtractionToSeparateImageFilter.h
Algorithms/itkTensorImageToDiffusionImageFilter.h
Algorithms/itkTensorToL2NormImageFilter.h
Algorithms/itkGaussianInterpolateImageFunction.h
Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h
Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h
Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h
Algorithms/itkCartesianToPolarVectorImageFilter.h
Algorithms/itkPolarToCartesianVectorImageFilter.h
Algorithms/itkDistanceMapFilter.h
Algorithms/itkProjectionFilter.h
Algorithms/itkResidualImageFilter.h
Algorithms/itkExtractChannelFromRgbaImageFilter.h
Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h
Algorithms/itkMergeDiffusionImagesFilter.h
Algorithms/itkDwiPhantomGenerationFilter.h
Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h
Algorithms/itkMrtrixPeakImageConverter.h
Algorithms/itkFslPeakImageConverter.h
Algorithms/itkShCoefficientImageImporter.h
Algorithms/itkShCoefficientImageExporter.h
Algorithms/itkOdfMaximaExtractionFilter.h
Algorithms/itkResampleDwiImageFilter.h
Algorithms/itkDwiGradientLengthCorrectionFilter.h
Algorithms/itkAdcImageFilter.h
Algorithms/itkDwiNormilzationFilter.h
Algorithms/itkSplitDWImageFilter.h
Algorithms/itkRemoveDwiChannelFilter.h
Algorithms/itkExtractDwiChannelFilter.h
Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h
- Algorithms/mitkDiffusionImageToDiffusionImageFilter.h
Algorithms/itkNonLocalMeansDenoisingFilter.h
Algorithms/itkVectorImageToImageFilter.h
)
set( TOOL_FILES
)