diff --git a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h index 4789f31764..06f9cee942 100644 --- a/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h +++ b/Modules/DiffusionImaging/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h @@ -1,145 +1,163 @@ /*=================================================================== 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 __mitkDiffusionImage__h #define __mitkDiffusionImage__h #include "mitkImage.h" #include "itkVectorImage.h" #include "itkVectorImageToImageAdaptor.h" +#include namespace mitk { - /** +/** * \brief this class encapsulates diffusion volumes (vectorimages not * yet supported by mitkImage) */ - template - class DiffusionImage : public Image - { +template +class DiffusionImage : public Image +{ - public: +public: typedef TPixelType PixelType; typedef typename itk::VectorImage - ImageType; + ImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, - GradientDirectionType > GradientDirectionContainerType; + GradientDirectionType > GradientDirectionContainerType; typedef itk::VectorImageToImageAdaptor< TPixelType, 3 > - AdaptorType; + AdaptorType; typedef vnl_matrix_fixed< double, 3, 3 > MeasurementFrameType; mitkClassMacro( DiffusionImage, Image ); itkNewMacro(Self); //void SetRequestedRegionToLargestPossibleRegion(); //bool RequestedRegionIsOutsideOfTheBufferedRegion(); //virtual bool VerifyRequestedRegion(); //void SetRequestedRegion(itk::DataObject *data); void AverageRedundantGradients(double precision); GradientDirectionContainerType::Pointer CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions); void CorrectDKFZBrokenGradientScheme(double precision); typename ImageType::Pointer GetVectorImage() { return m_VectorImage; } void SetVectorImage(typename ImageType::Pointer image ) { this->m_VectorImage = image; } void InitializeFromVectorImage(); void SetDisplayIndexForRendering(int displayIndex); GradientDirectionContainerType::Pointer GetDirections() { return m_Directions; } void SetDirections( GradientDirectionContainerType::Pointer directions ) { this->m_Directions = directions; } void SetDirections(const std::vector > directions) { - m_Directions = GradientDirectionContainerType::New(); - for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); - } + m_Directions = GradientDirectionContainerType::New(); + for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); + } } GradientDirectionContainerType::Pointer GetOriginalDirections() { return m_OriginalDirections; } void SetOriginalDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; this->ApplyMeasurementFrame(); } void SetOriginalDirections(const std::vector > directions) { - m_OriginalDirections = GradientDirectionContainerType::New(); - for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); - } - this->ApplyMeasurementFrame(); + m_OriginalDirections = GradientDirectionContainerType::New(); + for(unsigned int i=0; iInsertElement( i, directions[i].Get_vnl_vector() ); + } + this->ApplyMeasurementFrame(); } MeasurementFrameType GetMeasurementFrame() { return m_MeasurementFrame; } void SetMeasurementFrame( MeasurementFrameType mFrame ) { this->m_MeasurementFrame = mFrame; this->ApplyMeasurementFrame(); } itkGetMacro(B_Value, float); itkSetMacro(B_Value, float); float GetB_Value(int i) { - if(i > m_Directions->Size()-1) - return -1; - - if(m_Directions->ElementAt(i).one_norm() <= 0.0) - { - return 0; - } - else - { - double twonorm = m_Directions->ElementAt(i).two_norm(); - return m_B_Value*twonorm*twonorm ; - } + if(i > m_Directions->Size()-1) + return -1; + + if(m_Directions->ElementAt(i).one_norm() <= 0.0) + { + return 0; + } + else + { + double twonorm = m_Directions->ElementAt(i).two_norm(); + return m_B_Value*twonorm*twonorm ; + } } bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); int GetNumDirections(); int GetNumB0(); std::vector GetB0Indices(); bool IsMultiBval(); - protected: + void GetBvalueList(std::map > * bValMap) + { + if(bValMap != 0){ + GradientDirectionContainerType::ConstIterator gdcit; + for( gdcit = this->m_Directions->Begin(); gdcit != this->m_Directions->End(); ++gdcit) + { + float currentBvalue = std::floor(GetB_Value(gdcit.Index())); + double rounded = int((currentBvalue+7.5)/10)*10; + (*bValMap)[rounded].push_back(gdcit.Index()); + } + }else + { + MITK_ERROR << "Call method with 0 parameter"; + } + } + + +protected: DiffusionImage(); virtual ~DiffusionImage(); void ApplyMeasurementFrame(); typename ImageType::Pointer m_VectorImage; GradientDirectionContainerType::Pointer m_Directions; GradientDirectionContainerType::Pointer m_OriginalDirections; float m_B_Value; typename AdaptorType::Pointer m_VectorImageAdaptor; int m_DisplayIndex; MeasurementFrameType m_MeasurementFrame; - }; +}; } // namespace mitk #include "mitkDiffusionImage.txx" #endif /* __mitkDiffusionImage__h */ diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp new file mode 100644 index 0000000000..8bac58cbcf --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -0,0 +1,1135 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html 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 __itkDiffusionMultiShellQballReconstructionImageFilter_cpp +#define __itkDiffusionMultiShellQballReconstructionImageFilter_cpp + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +#include "mitkSphericalHarmonicsFunctions.h" +#include "mitkVNLVectorFunctions.h" +#include "itkPointShell.h" +#include + + +namespace itk { + +template< class T, class TG, class TO, int L, int NODF> +DiffusionMultiShellQballReconstructionImageFilter +::DiffusionMultiShellQballReconstructionImageFilter() : + m_GradientDirectionContainer(NULL), + m_NumberOfGradientDirections(0), + m_NumberOfBaselineImages(1), + m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()), + m_BValue(1.0), + m_Lambda(0.0), + m_IsHemisphericalArrangementOfGradientDirections(false), + m_IsArithmeticProgession(false), + m_ReconstructionType(Mode_Standard1Shell) +{ + // At least 1 inputs is necessary for a vector image. + // For images added one at a time we need at least six + this->SetNumberOfRequiredInputs( 1 ); +} + + +template +void DiffusionMultiShellQballReconstructionImageFilter +::Normalize( OdfPixelType & out) +{ + for(int i=0; i +double DiffusionMultiShellQballReconstructionImageFilter +::CalculateThreashold(const double value, const double delta) +{ + return (value<0)*(0.5*delta) + (value>=0 && value=delta && value<1-delta)*value+(value>=1-delta && value<1)*(1-0.5*delta-0.5*((1-value)*(1-value))/delta) + + (value>=1)*(1-0.5*delta); +} + +template +void DiffusionMultiShellQballReconstructionImageFilter +::Threshold(vnl_vector & vec, double delta) +{ + if (delta==0){ //Clip attenuation values. If att<0 => att=0, if att>1 => att=1 + for (int i=0; i=0 && vec[i]<=1)*vec[i]+(vec[i]>1); + } + else{ //Use function from Aganj et al, MRM, 2010 + for (int i=0; i< vec.size(); i++) + vec[i]=CalculateThreashold(vec[i], delta); + } +} + +template +void DiffusionMultiShellQballReconstructionImageFilter +::Threshold(vnl_matrix & mat, double delta) +{ + if (delta==0){ //Clip attenuation values. If att<0 => att=0, if att>1 => att=1 + for (int i=0; i=0 && mat(i,j)<=1)*mat(i,j)+(mat(i,j)>1); + } + else{ //Use function from Aganj et al, MRM, 2010 + for (int i=0; i +void DiffusionMultiShellQballReconstructionImageFilter +::Projection1( vnl_matrix & E, double delta ) +{ + + const double sF = sqrt(5.0); + + vnl_vector vOnes(E.rows()); + vOnes.fill(1.0); + + vnl_matrix T0(E.rows(), E.cols()); + vnl_matrix C(E.rows(), 7); + vnl_matrix A(E.rows(), 7); + vnl_matrix B(E.rows(), 7); + + vnl_vector s0(E.rows()); + vnl_vector a0(E.rows()); + vnl_vector b0(E.rows()); + vnl_vector ta(E.rows()); + vnl_vector tb(E.rows()); + vnl_vector e(E.rows()); + vnl_vector m(E.rows()); + vnl_vector a(E.rows()); + vnl_vector b(E.rows()); + + + // logarithmierung aller werte in E + for(int i = 0 ; i < E.rows(); i++) + { + for(int j = 0 ; j < E.cols(); j++) + { + T0(i,j) = -log(E(i,j)); + } + } + + + //T0 = -T0.apply(std::log); + + // Summeiere Zeilenweise über alle Shells sum = E1+E2+E3 + for(int i = 0 ; i < E.rows(); i++) + { + s0[i] = T0(i,0) + T0(i,1) + T0(i,2); + } + + + for(int i = 0; i < E.rows(); i ++) + { + // Alle Signal-Werte auf der Ersten shell E(N,0) normiert auf s0 + a0 = E(i,0) / s0[i]; + // Alle Signal-Werte auf der Zweiten shell E(N,1) normiert auf s0 + b0 = E(i,1) / s0[i]; + } + + ta = a0 * 3.0; + tb = b0 * 3.0; + e = tb - (ta * 2.0); + m = (tb * 2.0 ) + ta; + + for(int i = 0; i < E.rows(); i++) + { + C(i,0) = tb[i] < 1+3*delta && 0.5+1.5*(sF+1)*delta < ta[i] && ta[i] < 1-3* (sF+2) *delta; + C(i,1) = e[i] <= -1 +3*(2*sF+5)* delta && ta[i] >= 1-3*(sF+2)*delta; + C(i,2) = m[i] > 3 -3*sF*delta && -1+3*(2*sF+5)*delta < e[i] && e[i]<-3*sF*delta; + C(i,3) = m[i] >= 3-3*sF*delta && e[i] >= -3 *sF * delta; + C(i,4) = 2.5 + 1.5*(5+sF)*delta < m[i] && m[i] < 3-3*sF*delta && e[i] > -3*sF*delta; + C(i,5) = ta[i] <= 0.5+1.5 *(sF+1)*delta && m[i] <= 2.5 + 1.5 *(5+sF) * delta; + C(i,6) = !( C(i,0) || C(i,1) || C(i,2) || C(i,3) || C(i,4) || C(i,5) ); // ~ANY(C(i,[0-5] ),2) + + A(i,0)=(bool)C(i,0) * a0(i); + A(i,1)=(bool)C(i,1) * (1.0/3.0-(sF+2)*delta); + A(i,2)=(bool)C(i,2) * (0.2+0.8*a0(i)-0.4*b0(i)-delta/sF); + A(i,3)=(bool)C(i,3) * (0.2+delta/sF); + A(i,4)=(bool)C(i,4) * (0.2*a0(i)+0.4*b0(i)+2*delta/sF); + A(i,5)=(bool)C(i,5) * (1.0/6.0+0.5*(sF+1)*delta); + A(i,6)=(bool)C(i,6) * a0(i); + + B(i,0)=(bool)C(i,0) * (1.0/3.0+delta); + B(i,1)=(bool)C(i,1) * (1.0/3.0+delta); + B(i,2)=(bool)C(i,2) * (0.4-0.4*a0(i)+0.2*b0(i)-2*delta/sF); + B(i,3)=(bool)C(i,3) * (0.4-3*delta/sF); + B(i,4)=(bool)C(i,4) * (0.4*a0(i)+0.8*b0(i)-delta/sF); + B(i,5)=(bool)C(i,5) * (1.0/3.0+delta); + B(i,6)=(bool)C(i,6) * b0(i); + } + + + + for(int i = 0 ; i < E.rows(); i++) + { + double sumA = 0; + double sumB = 0; + for(int j = 0 ; j < 7; j++) + { + sumA += A(i,j); + sumB += B(i,j); + } + a[i] = sumA; + b[i] = sumB; + } + + for(int i = 0; i < E.rows(); i++) + { + E(i,0) = exp(-(a[i]*s0[i])); + E(i,1) = exp(-(b[i]*s0[i])); + E(i,2) = exp(-((1-a[i]-b[i])*s0[i])); + } + +} + + + + +template +void DiffusionMultiShellQballReconstructionImageFilter +::Projection2( vnl_vector & A, vnl_vector & a, vnl_vector & b, double delta0) +{ + + const double s6 = sqrt(6); + const double s15 = s6/2.0; + + vnl_vector delta(a.size()); + delta.fill(delta0); + + vnl_matrix AM(a.size(), 15); + vnl_matrix aM(a.size(), 15); + vnl_matrix bM(a.size(), 15); + vnl_matrix B(a.size(), 15); + + AM.set_column(0, A); + AM.set_column(1, A); + AM.set_column(2, A); + AM.set_column(3, delta); + AM.set_column(4, (A+a-b - (delta*s6))/3.0); + AM.set_column(5, delta); + AM.set_column(6, delta); + AM.set_column(7, delta); + AM.set_column(8, A); + AM.set_column(9, 0.2*(a*2+A-2*(s6+1)*delta)); + AM.set_column(10,0.2*(b*(-2)+A+2-2*(s6+1)*delta)); + AM.set_column(11, delta); + AM.set_column(12, delta); + AM.set_column(13, delta); + AM.set_column(14, 0.5-(1+s15)*delta); + + + aM.set_column(0, a); + aM.set_column(1, a); + aM.set_column(2, -delta + 1); + aM.set_column(3, a); + aM.set_column(4, (A*2+a*5+b+s6*delta)/6.0); + aM.set_column(5, a); + aM.set_column(6, -delta + 1); + aM.set_column(7, 0.5*(a+b)+(1+s15)*delta); + aM.set_column(8, -delta + 1); + aM.set_column(9, 0.2*(a*4+A*2+(s6+1)*delta)); + aM.set_column(10, -delta + 1); + aM.set_column(11, (s6+3)*delta); + aM.set_column(12, -delta + 1); + aM.set_column(13, -delta + 1); + aM.set_column(14, -delta + 1); + + bM.set_column(0, b); + bM.set_column(1, delta); + bM.set_column(2, b); + bM.set_column(3, b); + bM.set_column(4, (A*(-2)+a+b*5-s6*delta)/6.0); + bM.set_column(5, delta); + bM.set_column(6, b); + bM.set_column(7, 0.5*(a+b)-(1+s15)*delta); + bM.set_column(8, delta); + bM.set_column(9, delta); + bM.set_column(10, 0.2*(b*4-A*2+1-(s6+1)*delta)); + bM.set_column(11, delta); + bM.set_column(12, delta); + bM.set_column(13, -delta*(s6+3) + 1); + bM.set_column(14, delta); + + delta0 *= 0.99; + + for(int i = 0 ; i < a.size(); i ++) + { + for(int j = 0 ; j < 15; j ++) + { + B(i,j) = delta0 < AM(i,j) && 2 * (AM(i,j) + delta0 * s15) < aM(i,j) - bM(i,j) && bM(i,j) > delta0 && aM(i,j) < 1- delta0; + } + } + + vnl_matrix R2(a.size(), 15); + vnl_matrix A_(a.size(), 15); + vnl_matrix a_(a.size(), 15); + vnl_matrix b_(a.size(), 15); + + + + vnl_matrix OnesVecMat(1, 15); + OnesVecMat.fill(1.0); + + vnl_matrix AVecMat(a.size(), 1); + AVecMat.set_column(0,A); + + vnl_matrix aVecMat(a.size(), 1); + aVecMat.set_column(0,a); + + vnl_matrix bVecMat(a.size(), 1); + bVecMat.set_column(0,b); + + A_ = AM - (AVecMat * OnesVecMat); + a_ = aM - (aVecMat * OnesVecMat); + b_ = bM - (bVecMat * OnesVecMat); + + for(int i = 0 ; i < a.size(); i++) + for(int j = 0 ; j < 15; j++) + { + A_(i,j) *= A_(i,j); + a_(i,j) *= a_(i,j); + b_(i,j) *= b_(i,j); + } + + + R2 = A_ + a_ + b_; + + for(int i = 0 ; i < a.size(); i ++) + { + for(int j = 0 ; j < 15; j ++) + { + if(B(i,j) == 0) R2(i,j) = 1e20; + } + } + + std::vector indicies(a.size()); + + // suche den spalten-index der zu der kleinsten Zahl einer Zeile korrespondiert + for(int i = 0 ; i < a.size(); i++) + { + unsigned int index = 0; + double minvalue = 999; + for(int j = 0 ; j < 15 ; j++) + { + if(R2(i,j) < minvalue){ + minvalue = R2(i,j); + index = j; + } + } + indicies[i] = index; + } + + for(int i = 0 ; i < a.size(); i++) + { + A[i] = AM(i,indicies[i]); + a[i] = aM(i,indicies[i]); + b[i] = bM(i,indicies[i]); + } +} + + +template +void DiffusionMultiShellQballReconstructionImageFilter +::S_S0Normalization( vnl_vector & vec, typename NumericTraits::AccumulateType b0 ) +{ + + double b0f = (double)b0; + for(int i = 0; i < vec.size(); i++) + { + if (b0f==0) + b0f = 0.01; + if(vec[i] >= b0f) + vec[i] = b0f - 0.001; + vec[i] /= b0f; + } + +} + + +template +void DiffusionMultiShellQballReconstructionImageFilter +::S_S0Normalization( vnl_matrix & mat, typename NumericTraits::AccumulateType b0 ) +{ + double b0f = (double)b0; + for(int i = 0; i < mat.rows(); i++) + { + for( int j = 0; j < mat.cols(); j++ ){ + if (b0f==0) + b0f = 0.01; + if(mat(i,j) >= b0f) + mat(i,j) = b0f - 0.001; + mat(i,j) /= b0f; + } + } +} + + + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::DoubleLogarithm(vnl_vector & vec) +{ + for(int i = 0; i < vec.size(); i++) + { + vec[i] = log(-log(vec[i])); + } +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::SetGradientImage( GradientDirectionContainerType *gradientDirection + , const GradientImagesType *gradientImage + , float bvalue) +{ + this->m_BValue = bvalue; + this->m_GradientDirectionContainer = gradientDirection; + this->m_NumberOfBaselineImages = 0; + this->m_ReconstructionType = Mode_Standard1Shell; + + GradientDirectionContainerType::ConstIterator gdcit; + for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) + { + double bValueKey = int(((m_BValue* gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; + m_GradientIndexMap[bValueKey].push_back(gdcit.Index()); + } + + //if(listOfUserSelctedBValues.size() == 0){ + // itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : No list Of User Selcted B Values available"); + //} + if(m_GradientIndexMap.size() == 0){ + itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : no GradientIndexMapAvalible"); + } + //if(listOfUserSelctedBValues.size() != m_GradientIndexMap.size()){ + // itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : The number of user selected B Values != number of Image BValues"); + //} + + if(m_GradientIndexMap.size() == 4){ + + GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); + it++; + const int b1 = (*it).first; + it++; + const int b2 = (*it).first; + it++; + const int b3 = (*it).first; + + if(b2 - b1 == b1 && b3 - b2 == b1 ) + { + m_ReconstructionType = Mode_Analytical3Shells; + } + } + + if(m_GradientIndexMap.size() > 2 && m_ReconstructionType != Mode_Analytical3Shells) + { + m_ReconstructionType = Mode_NumericalNShells; + } + + this->m_NumberOfBaselineImages = m_GradientIndexMap[0].size(); + this->m_NumberOfGradientDirections = gradientDirection->Size() - this->m_NumberOfBaselineImages; + + // ensure that the gradient image we received has as many components as + // the number of gradient directions + if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections ) + { + itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages + << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages + << " directions specified but image has " << gradientImage->GetVectorLength() + << " components."); + } + + this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::BeforeThreadedGenerateData() +{ + itk::TimeProbe clock; + clock.Start(); + + if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) /* && m_NumberOfGradientDirections < 6 */ ) + { + itkExceptionMacro( << "At least " << ((L+1)*(L+2))/2 << " gradient directions are required" ); + } + // Input must be an itk::VectorImage. + std::string gradientImageClassName(this->ProcessObject::GetInput(0)->GetNameOfClass()); + + if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) + itkExceptionMacro( << "There is only one Gradient image. I expect that to be a VectorImage. But its of type: " << gradientImageClassName ); + + m_BZeroImage = BZeroImageType::New(); + typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing + m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin + m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction + m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); + m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); + m_BZeroImage->Allocate(); + + this->ComputeReconstructionMatrix(); + + clock.Stop(); + MITK_INFO << "Before GenerateData : " << clock.GetTotal(); + +} + + + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread) +{ + // Get output image pointer + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + + // ImageRegionIterator for the output image + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); + + // ImageRegionIterator for the BZero (output) image + ImageRegionConstIterator< BZeroImageType > bzeroImageIterator(m_BZeroImage, outputRegionForThread); + bzeroImageIterator.GoToBegin(); + + IndiciesVector SignalIndicies = m_GradientIndexMap[1]; + + // if the gradient directiosn aragement is hemispherical, duplicate all gradient directions + // alone, interested in the value, the direction can be neglected + if(m_IsHemisphericalArrangementOfGradientDirections){ + int NumbersOfGradientIndicies = SignalIndicies.size(); + for (int i = 0 ; i < NumbersOfGradientIndicies; i++) + SignalIndicies.push_back(SignalIndicies[i]); + } + + + // Get input gradient image pointer + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + + // Const ImageRegionIterator for input gradient image + typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; + GradientIteratorType git(gradientImagePointer, outputRegionForThread ); + git.GoToBegin(); + + typedef typename GradientImagesType::PixelType GradientVectorType; + + // iterate overall voxels of the gradient image region + while( ! git.IsAtEnd() ) + { + GradientVectorType b = git.Get(); + // ODF Vector + OdfPixelType odf(0.0); + + // Create the Signal Vector + vnl_vector SignalVector(m_NumberOfGradientDirections); + if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) + { + + for( unsigned int i = 0; i< SignalIndicies.size(); i++ ) + { + SignalVector[i] = static_cast(b[SignalIndicies[i]]); + } + + // apply threashold an generate ln(-ln(E)) signal + // Replace SignalVector with PreNormalized SignalVector + S_S0Normalization(SignalVector, bzeroImageIterator.Get()); + DoubleLogarithm(SignalVector); + + // ODF coeffs-vector + vnl_vector coeffs(m_NumberCoefficients); + + // approximate ODF coeffs + coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); + coeffs[0] = 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + + odf = mitk::mitk_vnl_function::element_cast(( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs )).data_block(); + odf *= (QBALL_ANAL_RECON_PI*4/NODF); + } + // set ODF to ODF-Image + oit.Set( odf ); + ++oit; + ++bzeroImageIterator; + ++git; + + } + + MITK_INFO << "One Thread finished reconstruction"; +} + +#include +//#include +//#include + + + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread) +{ + + + // vnl_levenberg_marquardt LMOptimizer = new vnl_levenberg_marquardt(); + +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread) +{ + typedef typename GradientImagesType::PixelType GradientVectorType; + + ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread); + bzeroIterator.GoToBegin(); + + IndiciesVector BZeroIndicies = m_GradientIndexMap[0]; + + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + + // Const ImageRegionIterator for input gradient image + typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; + GradientIteratorType git(gradientImagePointer, outputRegionForThread ); + git.GoToBegin(); + + + while( ! git.IsAtEnd() ) + { + GradientVectorType b = git.Get(); + + // compute the average bzero signal + typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; + for(unsigned int i = 0; i < BZeroIndicies.size(); ++i) + { + b0 += b[BZeroIndicies[i]]; + } + b0 /= BZeroIndicies.size(); + + bzeroIterator.Set(b0); + ++bzeroIterator; + ++git; + } +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread) +{ + + typedef typename GradientImagesType::PixelType GradientVectorType; + + // Input Gradient Image and Output ODF Image + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); + + + // Define Image iterators + ImageRegionIterator< OutputImageType > odfOutputImageIterator(outputImage, outputRegionForThread); + ImageRegionConstIterator< BZeroImageType > bzeroImageIterator(m_BZeroImage, outputRegionForThread); + ImageRegionConstIterator< GradientImagesType > gradientInputImageIterator(gradientImagePointer, outputRegionForThread ); + + // All iterators seht to Begin of the specific OutputRegion + odfOutputImageIterator.GoToBegin(); + bzeroImageIterator.GoToBegin(); + gradientInputImageIterator.GoToBegin(); + + // Get Shell Indicies for all non-BZero Gradients + // it MUST be a arithmetic progression eg.: 1000, 2000, 3000 + GradientIndexMapIteraotr it = m_GradientIndexMap.begin(); + it++; + // it = b-value = 1000 + IndiciesVector Shell1Indiecies = (*it).second; + it++; + // it = b-value = 2000 + IndiciesVector Shell2Indiecies = (*it).second; + it++; + // it = b-value = 3000 + IndiciesVector Shell3Indiecies = (*it).second; + + + // if input data is a hemispherical arragement, duplicate eache gradient for each shell + if(m_IsHemisphericalArrangementOfGradientDirections){ + int NumbersOfGradientIndicies = Shell1Indiecies.size(); + for (int i = 0 ; i < NumbersOfGradientIndicies; i++){ + Shell1Indiecies.push_back(Shell1Indiecies[i]); + Shell2Indiecies.push_back(Shell2Indiecies[i]); + Shell3Indiecies.push_back(Shell3Indiecies[i]); + } + } + + // Nx3 Signal Vector with E(0) = Shell 1, E(1) = Shell 2, E(2) = Shell 3 + vnl_matrix< double > * E = new vnl_matrix(Shell1Indiecies.size(), 3); + + vnl_vector * AlphaValues = new vnl_vector(Shell1Indiecies.size()); + vnl_vector * BetaValues = new vnl_vector(Shell1Indiecies.size()); + vnl_vector * LAValues = new vnl_vector(Shell1Indiecies.size()); + vnl_vector * PValues = new vnl_vector(Shell1Indiecies.size()); + + OdfPixelType odf(0.0); + + // iterate overall voxels of the gradient image region + while( ! gradientInputImageIterator.IsAtEnd() ) + { + if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) + { + // Get the Signal-Value for each Shell at each direction (specified in the ShellIndicies Vector .. this direction corresponse to this shell...) + GradientVectorType b = gradientInputImageIterator.Get(); + for(int i = 0 ; i < Shell1Indiecies.size(); i++) + { + E->put(i,0, static_cast(b[Shell1Indiecies[i]])); + E->put(i,1, static_cast(b[Shell2Indiecies[i]])); + E->put(i,2, static_cast(b[Shell3Indiecies[i]])); + } + + //Approximated-Signal by SH fit - using the specific shell directions and values + // approximated Signal : S = SHBasis * Coeffs + // with Coeffs : C = (B_T * B + lambda * L) ^ -1 * B_T * OS + // OS := Original-Signal + E->set_column(1, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(1)))); + E->set_column(2, (*m_SHBasisMatrix) * ((*m_SignalReonstructionMatrix) * (E->get_column(2)))); + + // Normalize the Signal: Si/S0 + S_S0Normalization(*E,bzeroImageIterator.Get()); + + //Implements Eq. [19] and Fig. 4. + Threshold(*E); + + //inqualities [31]. Taking the lograithm of th first tree inqualities + //convert the quadratic inqualities to linear ones. + Projection1(*E); + + double E1, E2, E3, P2,A,B2,B,P,alpha,beta,lambda, ER1, ER2; + + for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ ) + { + + E1 = E->get(i,0); + E2 = E->get(i,1); + E3 = E->get(i,2); + + P2 = E2-E1*E1; + A = (E3 -E1*E2) / ( 2* P2); + B2 = A * A -(E1 * E3 - E2 * E2) /P2; + B = 0; + if(B2 > 0) B = sqrt(B2); + + P = 0; + if(P2 > 0) P = sqrt(P2); + + alpha = A + B; + beta = A - B; + + + lambda = 0.5 + 0.5 * std::sqrt(1 - std::pow((P * 2 ) / (alpha - beta), 2));; + ER1 = std::fabs(lambda * (alpha - beta) + (beta - E1 )) + + std::fabs(lambda * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 )) + + std::fabs(lambda * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 )); + ER2 = std::fabs((lambda-1) * (alpha - beta) + (beta - E1 )) + + std::fabs((lambda-1) * (std::pow(alpha, 2) - std::pow(beta, 2)) + (std::pow(beta, 2) - E2 )) + + std::fabs((lambda-1) * (std::pow(alpha, 3) - std::pow(beta, 3)) + (std::pow(beta, 3) - E3 )); + + + PValues->put(i, P); + AlphaValues->put(i, alpha); + BetaValues->put(i, beta); + LAValues->put(i,(lambda * (ER1 < ER2)) + ((1-lambda) * (ER1 >= ER2))); + + } + Projection2(*PValues, *AlphaValues, *BetaValues); + + //Threshold(*AlphaValues); + //Threshold(*BetaValues); + + DoubleLogarithm(*AlphaValues); + DoubleLogarithm(*BetaValues); + + vnl_vector SignalVector(element_product((*LAValues) , (*AlphaValues)-(*BetaValues)) + (*BetaValues)); + vnl_vector coeffs((*m_CoeffReconstructionMatrix) *SignalVector ); + + // the first coeff is a fix value + coeffs[0] = 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); + + // Cast the Signal-Type from double to float for the ODF-Image + odf = mitk::mitk_vnl_function::element_cast( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); + odf *= (QBALL_ANAL_RECON_PI*4/NODF); + //Normalize(odf); + } + + // set ODF to ODF-Image + odfOutputImageIterator.Set( odf ); + ++odfOutputImageIterator; + // iterate + ++bzeroImageIterator; + ++gradientInputImageIterator; + } + + + + MITK_INFO << "THREAD FINISHED"; + delete E; + delete AlphaValues; + delete BetaValues; + delete PValues; + delete LAValues; +} + + +template< class T, class TG, class TO, int L, int NODF> +vnl_vector DiffusionMultiShellQballReconstructionImageFilter +::AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1Indicies,const IndiciesVector * shell2Indicies ,const IndiciesVector * shell3Indicies, vnl_vector) +{ + + + +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads) +{ + itk::TimeProbe clock; + + GenerateAveragedBZeroImage(outputRegionForThread); + + clock.Start(); + switch(m_ReconstructionType) + { + case Mode_Standard1Shell: + StandardOneShellReconstruction(outputRegionForThread); + break; + case Mode_Analytical3Shells: + AnalyticalThreeShellReconstruction(outputRegionForThread); + break; + case Mode_NumericalNShells: + break; + } + clock.Stop(); + MITK_INFO << "Reconstruction in : " << clock.GetTotal() << " TU"; + + + +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter:: +ComputeSphericalHarmonicsBasis(vnl_matrix * QBallReference, vnl_matrix *SHBasisOutput, vnl_matrix* LaplaciaBaltramiOutput, vnl_vector* SHOrderAssociation, vnl_matrix* SHEigenvalues ) +{ + for(unsigned int i=0; i< (*SHBasisOutput).rows(); i++) + { + for(int k = 0; k <= L; k += 2) + { + for(int m =- k; m <= k; m++) + { + int j = ( k * k + k + 2 ) / 2 + m - 1; + + // Compute SHBasisFunctions + double phi = (*QBallReference)(0,i); + double th = (*QBallReference)(1,i); + (*SHBasisOutput)(i,j) = mitk::sh::Yj(m,k,th,phi); + + // Laplacian Baltrami Order Association + if(LaplaciaBaltramiOutput) + (*LaplaciaBaltramiOutput)(j,j) = k*k*(k + 1)*(k+1); + + // SHEigenvalues with order Accosiation kj + if(SHEigenvalues) + (*SHEigenvalues)(j,j) = -k* (k+1); + + // Order Association + if(SHOrderAssociation) + (*SHOrderAssociation)[j] = k; + + } + } + } +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::ComputeFunkRadonTransformationMatrix(vnl_vector* SHOrderAssociationReference, vnl_matrix* FRTMatrixOutput ) +{ + for(int i=0; i +bool DiffusionMultiShellQballReconstructionImageFilter +::CheckHemisphericalArrangementOfGradientDirections() +{ + // handle acquisition schemes where only half of the spherical + // shell is sampled by the gradient directions. In this case, + // each gradient direction is duplicated in negative direction. + vnl_vector centerMass(3); + centerMass.fill(0.0); + int count = 0; + GradientDirectionContainerType::ConstIterator gdcit1; + for( gdcit1 = this->m_GradientDirectionContainer->Begin(); gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1) + { + if(gdcit1.Value().one_norm() > 0.0) + { + centerMass += gdcit1.Value(); + count ++; + } + } + centerMass /= count; + if(centerMass.two_norm() > 0.1) + { + return false; + } + return true; +} + +template< class T, class TG, class TO, int L, int NOdfDirections> +void DiffusionMultiShellQballReconstructionImageFilter +::ComputeReconstructionMatrix() +{ + + typedef std::auto_ptr< vnl_matrix< double> > MatrixDoublePtr; + typedef std::auto_ptr< vnl_vector< int > > VectorIntPtr; + typedef std::auto_ptr< vnl_matrix_inverse< double > > InverseMatrixDoublePtr; + + std::map >::const_iterator it = (m_GradientIndexMap.begin()); + it++; + const std::vector gradientIndiciesVector= (*it).second; + + int numberOfGradientDirections = gradientIndiciesVector.size(); + + + if( numberOfGradientDirections < (((L+1)*(L+2))/2) || numberOfGradientDirections < 6 ) + { + itkExceptionMacro( << "At least (L+1)(L+2)/2 gradient directions for each shell are required; current : " << numberOfGradientDirections ); + } + + + CheckDuplicateDiffusionGradients(); + + // check if gradient directions are arrangement as a hemisphere(true) or sphere(false) + m_IsHemisphericalArrangementOfGradientDirections = CheckHemisphericalArrangementOfGradientDirections(); + if(m_IsHemisphericalArrangementOfGradientDirections) numberOfGradientDirections *= 2; + + + MatrixDoublePtr Q(new vnl_matrix(3, numberOfGradientDirections)); + Q->fill(0.0); + // Cartesian to spherical coordinates + { + int j = 0; + + for(int i = 0; i < gradientIndiciesVector.size(); i++) + { + + double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); + double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); + double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); + double cart[3]; + mitk::sh::Cart2Sph(x,y,z,cart); + (*Q)(0,j) = cart[0]; + (*Q)(1,j) = cart[1]; + (*Q)(2,j++) = cart[2]; + + } + if(m_IsHemisphericalArrangementOfGradientDirections) + { + for(int i = 0; i < gradientIndiciesVector.size(); i++) + { + double x = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(0); + double y = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(1); + double z = this->m_GradientDirectionContainer->ElementAt(gradientIndiciesVector[i]).get(2); + double cart[3]; + mitk::sh::Cart2Sph(x,y,z,cart); + (*Q)(0,j) = cart[0]; + (*Q)(1,j) = cart[1]; + (*Q)(2,j++) = cart[2]; + } + } + } + + const int LOrder = L; + m_NumberCoefficients = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; + MITK_INFO << m_NumberCoefficients; + m_SHBasisMatrix = new vnl_matrix(numberOfGradientDirections,m_NumberCoefficients); + m_SHBasisMatrix->fill(0.0); + VectorIntPtr SHOrderAssociation(new vnl_vector(m_NumberCoefficients)); + SHOrderAssociation->fill(0.0); + MatrixDoublePtr LaplacianBaltrami(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + LaplacianBaltrami->fill(0.0); + MatrixDoublePtr FRTMatrix(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + FRTMatrix->fill(0.0); + MatrixDoublePtr SHEigenvalues(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + SHEigenvalues->fill(0.0); + + + + // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector + ComputeSphericalHarmonicsBasis(Q.get() ,m_SHBasisMatrix, LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); + + // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj + ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get()); + + MatrixDoublePtr temp(new vnl_matrix(((m_SHBasisMatrix->transpose()) * (*m_SHBasisMatrix)) + (m_Lambda * (*LaplacianBaltrami)))); + + InverseMatrixDoublePtr pseudo_inv(new vnl_matrix_inverse((*temp))); + MatrixDoublePtr inverse(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + inverse->fill(0.0); + (*inverse) = pseudo_inv->inverse(); + + // ODF Factor ( missing 1/4PI ?? ) + double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); + + m_SignalReonstructionMatrix = new vnl_matrix((*inverse) * (m_SHBasisMatrix->transpose())); + + m_CoeffReconstructionMatrix = new vnl_matrix(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) )); + + // this code goes to the image adapter coeffs->odfs later + + vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); + + m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients); + m_ODFSphericalHarmonicBasisMatrix->fill(0.0); + + for(int i=0; i +template< class VNLType > +void DiffusionMultiShellQballReconstructionImageFilter +::printMatrix( VNLType * mat ) +{ + + std::stringstream stream; + + for(int i = 0 ; i < mat->rows(); i++) + { + stream.str(""); + for(int j = 0; j < mat->cols(); j++) + { + stream << (*mat)(i,j) << " "; + } + + } + MITK_INFO << stream.str(); +} + +template< class T, class TG, class TO, int L, int NODF> +bool DiffusionMultiShellQballReconstructionImageFilter +::CheckDuplicateDiffusionGradients() +{ + bool value = false; + + GradientIndexMapIteraotr mapIterator = m_GradientIndexMap.begin(); + while(mapIterator != m_GradientIndexMap.end()) + { + std::vector::const_iterator it1 = mapIterator->second.begin(); + std::vector::const_iterator it2 = mapIterator->second.begin(); + + for(; it1 != mapIterator->second.end(); ++it1) + { + for(; it2 != mapIterator->second.end(); ++it2) + { + if(m_GradientDirectionContainer->ElementAt(*it1) == m_GradientDirectionContainer->ElementAt(*it2) && it1 != it2) + { + itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." ); + value = true; + } + } + } + ++mapIterator; + } + return value; +} + + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::PrintSelf(std::ostream& os, Indent indent) const +{ + std::locale C("C"); + std::locale originalLocale = os.getloc(); + os.imbue(C); + + Superclass::PrintSelf(os,indent); + + //os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl; + if ( m_GradientDirectionContainer ) + { + os << indent << "GradientDirectionContainer: " + << m_GradientDirectionContainer << std::endl; + } + else + { + os << indent << + "GradientDirectionContainer: (Gradient directions not set)" << std::endl; + } + os << indent << "NumberOfGradientDirections: " << + m_NumberOfGradientDirections << std::endl; + os << indent << "NumberOfBaselineImages: " << + m_NumberOfBaselineImages << std::endl; + os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl; + os << indent << "BValue: " << m_BValue << std::endl; + + os.imbue( originalLocale ); +} + + + + +} + +#endif // __itkDiffusionMultiShellQballReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h new file mode 100644 index 0000000000..3e381c8221 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -0,0 +1,216 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html 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 __itkDiffusionMultiShellQballReconstructionImageFilter_h_ +#define __itkDiffusionMultiShellQballReconstructionImageFilter_h_ + +#include "itkImageToImageFilter.h" +#include "vnl/vnl_vector_fixed.h" +#include "vnl/vnl_matrix.h" +#include "vnl/algo/vnl_svd.h" +#include "itkVectorContainer.h" +#include "itkVectorImage.h" +#include + + + +namespace itk{ +/** \class DiffusionMultiShellQballReconstructionImageFilter + Aganj_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 Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; + typedef TReferenceImagePixelType ReferencePixelType; + typedef TGradientImagePixelType GradientPixelType; + typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; + typedef typename Superclass::InputImageType ReferenceImageType; + typedef Image< OdfPixelType, 3 > OdfImageType; + typedef OdfImageType OutputImageType; + typedef TOdfPixelType BZeroPixelType; + typedef Image< BZeroPixelType, 3 > BZeroImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + /** Typedef defining one (of the many) gradient images. */ + typedef Image< GradientPixelType, 3 > GradientImageType; + /** An alternative typedef defining one (of the many) gradient images. + * It will be assumed that the vectorImage has the same dimension as the + * Reference image and a vector length parameter of \c n (number of + * gradient directions)*/ + typedef VectorImage< GradientPixelType, 3 > GradientImagesType; + /** Holds the ODF reconstruction matrix */ + typedef vnl_matrix< TOdfPixelType >* OdfReconstructionMatrixType; + typedef vnl_matrix< double > * CoefficientMatrixType; + /** Holds each magnetic field gradient used to acquire one DWImage */ + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + /** Container to hold gradient directions of the 'n' DW measurements */ + typedef VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; + + typedef std::map > GradientIndexMap; + typedef std::map >::iterator GradientIndexMapIteraotr; + typedef std::vector IndiciesVector; + + // --------------------------------------------------------------------------------------------// + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter); + + /** set method to add gradient directions and its corresponding + * image. The image here is a VectorImage. The user is expected to pass the + * gradient directions in a container. The ith element of the container + * corresponds to the gradient direction of the ith component image the + * VectorImage. For the baseline image, a vector of all zeros + * should be set.*/ + void SetGradientImage( GradientDirectionContainerType *, const GradientImagesType *image , float bvalue);//, std::vector listOfUserSelctedBValues ); + + + /** Get reference image */ + virtual ReferenceImageType * GetReferenceImage() + { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); } + + /** Return the gradient direction. idx is 0 based */ + virtual GradientDirectionType GetGradientDirection( unsigned int idx) const + { + if( idx >= m_GradientDirectionContainer->Size() ) + { + itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); + } + return m_GradientDirectionContainer->ElementAt( idx+1 ); + } + + void Normalize(OdfPixelType & odf ); + + void S_S0Normalization( vnl_vector & vec, typename NumericTraits::AccumulateType b0 = 0 ); + void S_S0Normalization( vnl_matrix & mat, typename NumericTraits::AccumulateType b0 = 0 ); + + void DoubleLogarithm(vnl_vector & vec); + + void Threshold(vnl_vector & vec, double delta = 0.01); + void Threshold(vnl_matrix & mat, double delta = 0.01); + double CalculateThreashold(const double value, const double delta); + + void Projection1( vnl_matrix & mat, double delta = 0.01); + void Projection2( vnl_vector & A, vnl_vector & alpha, vnl_vector & beta, double delta = 0.01); + + /** 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( BZeroImage, typename BZeroImageType::Pointer); + //itkGetMacro( ODFSumImage, typename BlaImage::Pointer); + + itkSetMacro( Lambda, double ); + itkGetMacro( Lambda, double ); + + itkGetConstReferenceMacro( BValue, TOdfPixelType); + + +protected: + DiffusionMultiShellQballReconstructionImageFilter(); + ~DiffusionMultiShellQballReconstructionImageFilter() {}; + void PrintSelf(std::ostream& os, Indent indent) const; + + void ComputeReconstructionMatrix(); + bool CheckDuplicateDiffusionGradients(); + void ComputeSphericalHarmonicsBasis(vnl_matrix* QBallReference, vnl_matrix* SHBasisOutput, vnl_matrix* LaplaciaBaltramiOutput, vnl_vector* SHOrderAssociation , vnl_matrix * SHEigenvalues); + void ComputeFunkRadonTransformationMatrix(vnl_vector* SHOrderAssociationReference, vnl_matrix* FRTMatrixOutput ); + bool CheckHemisphericalArrangementOfGradientDirections(); + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads ); + + vnl_vector AnalyticalThreeShellParameterEstimation(const IndiciesVector * shell1, const IndiciesVector * shell2, const IndiciesVector * shell3, vnl_vector b); + void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread); + void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread); + void NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread); + void GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread); + +private: + + enum ReconstructionType + { + Mode_Analytical3Shells, + Mode_NumericalNShells, + Mode_Standard1Shell + }; + + + + //CoefficientMatrixType m_ReconstructionMatrix; + CoefficientMatrixType m_CoeffReconstructionMatrix; + CoefficientMatrixType m_ODFSphericalHarmonicBasisMatrix; + CoefficientMatrixType m_SignalReonstructionMatrix; + CoefficientMatrixType m_SHBasisMatrix; + + /** 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; + + /** LeBihan's b-value for normalizing tensors */ + float m_BValue; + + typename BZeroImageType::Pointer m_BZeroImage; + + GradientIndexMap m_GradientIndexMap; + + double m_Lambda; + + bool m_IsHemisphericalArrangementOfGradientDirections; + + bool m_IsArithmeticProgession; + + int m_NumberCoefficients; + + ReconstructionType m_ReconstructionType; + + + + template< class VNLType > + void printMatrix( VNLType * mat ); + + + +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" +#endif + +#endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_ + diff --git a/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp new file mode 100644 index 0000000000..3f931d4374 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/mitkSphericalHarmonicsFunctions.cpp @@ -0,0 +1,91 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html 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. + +=========================================================================*/ + +#include + +#define _USE_MATH_DEFINES +#include + +#if BOOST_VERSION / 100000 > 0 +#if BOOST_VERSION / 100 % 1000 > 34 +#include +#include +#endif +#endif + +#include "mitkVector.h" +#include "mitkSphericalHarmonicsFunctions.h" +namespace mitk{ +using namespace ::boost::math; + +double mitk::sh::factorial(int number) { + if(number <= 1) return 1; + double result = 1.0; + for(int i=1; i<=number; i++) + result *= i; + return result; +} + +void mitk::sh::Cart2Sph(double x, double y, double z, double *cart) +{ + double phi, th, rad; + rad = sqrt(x*x+y*y+z*z); + if( rad < mitk::eps ) + { + th = M_PI/2; + phi = M_PI/2; + } + else + { + th = acos(z/rad); + phi = atan2(y, x); + } + cart[0] = phi; + cart[1] = th; + cart[2] = rad; +} + +double mitk::sh::legendre0(int l) +{ + if( l%2 != 0 ) + { + return 0; + } + else + { + double prod1 = 1.0; + for(int i=1;i + +#define _USE_MATH_DEFINES +#include + +#if BOOST_VERSION / 100000 > 0 +#if BOOST_VERSION / 100 % 1000 > 34 +#include +#endif +#endif + +#include "mitkSphericalHarmonicsFunctions.h" +namespace mitk{ + +#define SPHERICAL_HARMONICS_PI M_PI + +double mitk::ShericalHarmonicsFunctions::factorial(int number) { + if(number <= 1) return 1; + double result = 1.0; + for(int i=1; i<=number; i++) + result *= i; + return result; +} + +void mitk::ShericalHarmonicsFunctions::Cart2Sph(double x, double y, double z, double *cart) +{ + double phi, th, rad; + rad = sqrt(x*x+y*y+z*z); + th = atan2(z,sqrt(x*x+y*y)); + phi = atan2(y,x); + th = -th+SPHERICAL_HARMONICS_PI/2; + phi = -phi+SPHERICAL_HARMONICS_PI; + cart[0] = phi; + cart[1] = th; + cart[2] = rad; +} + +double mitk::ShericalHarmonicsFunctions::legendre0(int l) +{ + if( l%2 != 0 ) + { + return 0; + } + else + { + double prod1 = 1.0; + for(int i=1;i SPHERICAL_HARMONICS_PI ) + { + std::cout << "bad range" << std::endl; + return 0; + } + + if( phi < 0 || phi > 2*SPHERICAL_HARMONICS_PI ) + { + std::cout << "bad range" << std::endl; + return 0; + } + + double pml = 0; + double fac1 = factorial(l+m); + double fac2 = factorial(l-m); + + if( m<0 ) + { +#if BOOST_VERSION / 100000 > 0 +#if BOOST_VERSION / 100 % 1000 > 34 + pml = ::boost::math::legendre_p(l, -m, cos(theta)); +#else + std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; +#endif +#else + std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; +#endif + double mypow = pow(-1.0,-m); + double myfac = (fac1/fac2); + pml *= mypow*myfac; + } + else + { +#if BOOST_VERSION / 100000 > 0 +#if BOOST_VERSION / 100 % 1000 > 34 + pml = ::boost::math::legendre_p(l, m, cos(theta)); +#endif +#endif + } + + //std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl; + + double retval = sqrt(((2.0*(double)l+1.0)/(4.0*SPHERICAL_HARMONICS_PI))*(fac2/fac1)) * pml; + if( !complexPart ) + { + retval *= cos(m*phi); + } + else + { + retval *= sin(m*phi); + } + //std::cout << retval << std::endl; + return retval; +} + +double mitk::ShericalHarmonicsFunctions::Yj(int m, int k, double theta, double phi) +{ + if( -k <= m && m < 0 ) + { + return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false); + } + + if( m == 0 ) + return spherical_harmonic(0,k,theta,phi,false); + + if( 0 < m && m <= k ) + { + return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true); + } + + return 0; +} +} diff --git a/Modules/DiffusionImaging/Reconstruction/mitkVNLVectorFunctions.h b/Modules/DiffusionImaging/Reconstruction/mitkVNLVectorFunctions.h new file mode 100644 index 0000000000..15149936a2 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/mitkVNLVectorFunctions.h @@ -0,0 +1,145 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $ +Version: $Revision: 18127 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html 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 __mitkVNLVectorFunctions_h_ +#define __mitkVNLVectorFunctions_h_ + +#include "MitkDiffusionImagingExports.h" +#include "vnl/vnl_vector.h" + +namespace mitk +{ + +namespace mitk_vnl_function +{ + +template +vnl_vector element_abs (vnl_vector const& v1) +{ + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + result[i] = std::fabs(v1[i]); + + return result; +} + +template +vnl_vector element_cast (vnl_vector const& v1) +{ + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + result[i] = static_cast(v1[i]); + + return result; +} + + +template +vnl_vector element_pow (vnl_vector const& v1, int p) +{ + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + result[i] = std::pow(v1[i], p); + + return result; +} + +template +vnl_vector element_condition_smallerThan (vnl_vector const& v1, vnl_vector const& v2, vnl_vector const& v3) +{ + +#ifndef NDEBUG + if (v1.size() != v2.size()) + vnl_error_vector_dimension ("element_condition_smaller: condition value", v1.size(), v2.size()); + if (v1.size() != v3.size()) + vnl_error_vector_dimension ("element_condition_smaller: accept value", v1.size(), v3.size()); +#endif + + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + v1[i] < v2[i] ? result[i] = v3[i] : result[i] = 0; + + return result; +} + +template +vnl_vector element_condition_smallerThanEqual (vnl_vector const& v1, vnl_vector const& v2, vnl_vector const& v3) +{ + +#ifndef NDEBUG + if (v1.size() != v2.size()) + vnl_error_vector_dimension ("element_condition_smallerThanEqual: condition value", v1.size(), v2.size()); + if (v1.size() != v3.size()) + vnl_error_vector_dimension ("element_condition_smallerThanEqual: accept value", v1.size(), v3.size()); +#endif + + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + v1[i] <= v2[i] ? result[i] = v3[i] : result[i] = 0; + + return result; +} + +template +vnl_vector element_condition_greaterThan (vnl_vector const& v1, vnl_vector const& v2, vnl_vector const& v3) +{ + +#ifndef NDEBUG + if (v1.size() != v2.size()) + vnl_error_vector_dimension ("element_condition_greaterThan: condition value", v1.size(), v2.size()); + if (v1.size() != v3.size()) + vnl_error_vector_dimension ("element_condition_greaterThan: accept value", v1.size(), v3.size()); +#endif + + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + v1[i] > v2[i] ? result[i] = v3[i] : result[i] = 0; + + return result; +} + +template +vnl_vector element_condition_greaterThanEqual (vnl_vector const& v1, vnl_vector const& v2, vnl_vector const& v3) +{ + +#ifndef NDEBUG + if (v1.size() != v2.size()) + vnl_error_vector_dimension ("element_condition_greaterThanEqual: condition value", v1.size(), v2.size()); + if (v1.size() != v3.size()) + vnl_error_vector_dimension ("element_condition_greaterThanEqual: accept value", v1.size(), v3.size()); +#endif + + vnl_vector result(v1.size()); + + for(int i = 0 ; i < v1.size(); i++) + v1[i] >= v2[i] ? result[i] = v3[i] : result[i] = 0; + + return result; +} + + +} +} + +#endif //__mitkVNLVectorFunctions_h_ + diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index eac860f869..efddb2e1f3 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,219 +1,226 @@ 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 # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # DataStructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp Rendering/mitkConnectomicsNetworkMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/mitkTractAnalyzer.cpp # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.cpp Algorithms/Connectomics/mitkConnectomicsHistogramBase.cpp Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.cpp Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.cpp Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.cpp Algorithms/Connectomics/mitkConnectomicsHistogramCache.cpp Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.cpp # Tractography Tractography/itkStochasticTractographyFilter.h + + # Reconstruction + Reconstruction/mitkSphericalHarmonicsFunctions.cpp + ) set(H_FILES # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h Rendering/mitkConnectomicsNetworkMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h + Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h + Reconstruction/mitkSphericalHarmonicsFunctions.h + Reconstruction/mitkVNLVectorFunctions.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/TbssImages/mitkTbssImporter.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Datastructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.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/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h Algorithms/itkReduceDirectionGradientsFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.h Algorithms/Connectomics/mitkConnectomicsHistogramBase.h Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.h Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.h Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.h Algorithms/Connectomics/mitkConnectomicsHistogramCache.h Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h ) set( TOOL_FILES ) if(WIN32) endif(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )