diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp index 7bd21ff93e..d6d89badbe 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.cpp @@ -1,1128 +1,1066 @@ /*=================================================================== 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_cpp #define __itkDiffusionMultiShellQballReconstructionImageFilter_cpp #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionFunctionCollection.h" #include "itkPointShell.h" #include +#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) + m_ReconstructionType(Mode_Standard1Shell), + m_Interpolation_Flag(false), + m_Interpolation_SHT1_inv(0), + m_Interpolation_SHT2_inv(0), + m_Interpolation_SHT3_inv(0), + m_Interpolation_TARGET_SH(0) { // 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) +::Projection1(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) +double DiffusionMultiShellQballReconstructionImageFilter +::CalculateThreashold(const double value, const 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=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 -::Projection1( vnl_matrix & E, double delta ) +::Projection2( vnl_vector & E1,vnl_vector & E2, vnl_vector & E3, double delta ) { const double sF = sqrt(5.0); - vnl_vector vOnes(E.rows()); + vnl_vector vOnes(m_MaxDirections); 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_matrix T0(m_MaxDirections, 3); + vnl_matrix C(m_MaxDirections, 7); + vnl_matrix A(m_MaxDirections, 7); + vnl_matrix B(m_MaxDirections, 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()); + vnl_vector s0(m_MaxDirections); + vnl_vector a0(m_MaxDirections); + vnl_vector b0(m_MaxDirections); + vnl_vector ta(m_MaxDirections); + vnl_vector tb(m_MaxDirections); + vnl_vector e(m_MaxDirections); + vnl_vector m(m_MaxDirections); + vnl_vector a(m_MaxDirections); + vnl_vector b(m_MaxDirections); // logarithmierung aller werte in E - for(int i = 0 ; i < E.rows(); i++) + for(int i = 0 ; i < m_MaxDirections; i++) { - for(int j = 0 ; j < E.cols(); j++) - { - T0(i,j) = -log(E(i,j)); - } + T0(i,0) = -log(E1(i)); + T0(i,1) = -log(E2(i)); + T0(i,2) = -log(E3(i)); } - //T0 = -T0.apply(std::log); + // Summeiere Zeilenweise über alle Shells sum = E1+E2+E3 - for(int i = 0 ; i < E.rows(); i++) + for(int i = 0 ; i < m_MaxDirections; i++) { s0[i] = T0(i,0) + T0(i,1) + T0(i,2); } - - for(int i = 0; i < E.rows(); i ++) + for(int i = 0; i < m_MaxDirections; i ++) { // Alle Signal-Werte auf der Ersten shell E(N,0) normiert auf s0 - a0 = E(i,0) / s0[i]; + a0[i] = T0(i,0) / s0[i]; // Alle Signal-Werte auf der Zweiten shell E(N,1) normiert auf s0 - b0 = E(i,1) / s0[i]; + b0[i] = T0(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++) + for(int i = 0; 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) + 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= 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) = !((bool) C(i,0) ||(bool) C(i,1) ||(bool) C(i,2) ||(bool) C(i,3) ||(bool) C(i,4) ||(bool) 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++) + for(int i = 0 ; i < m_MaxDirections; 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++) + for(int i = 0; i < m_MaxDirections; 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])); + E1(i) = exp(-(a[i]*s0[i])); + E2(i) = exp(-(b[i]*s0[i])); + E3(i) = exp(-((1-a[i]-b[i])*s0[i])); } } - - - template void DiffusionMultiShellQballReconstructionImageFilter -::Projection2( vnl_vector & A, vnl_vector & a, vnl_vector & b, double delta0) +::Projection3( vnl_vector & A, vnl_vector & a, vnl_vector & b, double delta0) { const double s6 = sqrt(6.0); 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); + std::vector I(a.size()); + + for (int i=0; idelta0 && aM(i,j)<1-delta0) + R2(i,j) = (AM(i,j)-A(i))*(AM(i,j)-A(i))+ (aM(i,j)-a(i))*(aM(i,j)-a(i))+(bM(i,j)-b(i))*(bM(i,j)-b(i)); + else + R2(i,j) = 1e20; } - - - 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; + I[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]); + for (int i=0; i < A.size(); i++){ + A(i) = AM(i,(int)I[i]); + a(i) = aM(i,(int)I[i]); + b(i) = bM(i,(int)I[i]); } -} +} template void DiffusionMultiShellQballReconstructionImageFilter -::S_S0Normalization( vnl_vector & vec, typename NumericTraits::AccumulateType b0 ) +::S_S0Normalization( vnl_vector & vec, double S0 ) { - - 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; + if (S0==0) + S0 = 0.01; + vec[i] /= S0; } } - -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; + m_BValue = bvalue; + m_GradientDirectionContainer = gradientDirection; + m_NumberOfBaselineImages = 0; + if(m_BValueMap.size() == 0){ itkWarningMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : no GradientIndexMapAvalible"); GradientDirectionContainerType::ConstIterator gdcit; - for( gdcit = this->m_GradientDirectionContainer->Begin(); gdcit != this->m_GradientDirectionContainer->End(); ++gdcit) + for( gdcit = m_GradientDirectionContainer->Begin(); gdcit != m_GradientDirectionContainer->End(); ++gdcit) { double bValueKey = int(((m_BValue * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; m_BValueMap[bValueKey].push_back(gdcit.Index()); } } - if(m_BValueMap.find(0) == m_BValueMap.end()) { itkExceptionMacro(<< "DiffusionMultiShellQballReconstructionImageFilter.cpp : GradientIndxMap with no b-Zero indecies found: check input image"); } + + m_NumberOfBaselineImages = m_BValueMap[0].size(); + m_NumberOfGradientDirections = gradientDirection->Size() - m_NumberOfBaselineImages; + // ensure that the gradient image we received has as many components as + // the number of gradient directions + if( gradientImage->GetVectorLength() != m_NumberOfBaselineImages + m_NumberOfGradientDirections ) + { + itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << m_NumberOfBaselineImages + << "baselines = " << m_NumberOfGradientDirections + m_NumberOfBaselineImages + << " directions specified but image has " << gradientImage->GetVectorLength() + << " components."); + } + + + ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); + + std::string gradientImageClassName(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 * >( 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(); + +} + +template< class T, class TG, class TO, int L, int NODF> +void DiffusionMultiShellQballReconstructionImageFilter +::BeforeThreadedGenerateData() +{ + m_ReconstructionType = Mode_Standard1Shell; + if(m_BValueMap.size() == 4 ){ BValueMapIteraotr it = m_BValueMap.begin(); + it++; // skip b0 entry + const int b1 = it->first; + const int vecSize1 = it->second.size(); + IndiciesVector shell1 = it->second; it++; - const int b1 = (*it).first; - it++; - const int b2 = (*it).first; + const int b2 = it->first; + const int vecSize2 = it->second.size(); + IndiciesVector shell2 = it->second; it++; - const int b3 = (*it).first; + const int b3 = it->first; + const int vecSize3 = it->second.size(); + IndiciesVector shell3 = it->second; + // arithmetic progrssion if(b2 - b1 == b1 && b3 - b2 == b1 ) { + // check if Interpolation is needed + // if shells with different numbers of directions exist + m_Interpolation_Flag = false; + if(vecSize1 != vecSize2 || vecSize2 != vecSize3 || vecSize1 != vecSize3) + { + m_Interpolation_Flag = true; + MITK_INFO << "Shell interpolation: shells with different numbers of directions"; + }else // if each shell holds same numbers of directions, but the gradient direction differ more than one 1 degree + { + m_Interpolation_Flag = CheckForDifferingShellDirections(); + if(m_Interpolation_Flag) MITK_INFO << "Shell interpolation: gradient direction differ more than one 1 degree"; + } + m_ReconstructionType = Mode_Analytical3Shells; + + if(m_Interpolation_Flag) + { + IndiciesVector min_shell; + IndiciesVector max_shell; + int Interpolation_SHOrder = 10; + + //fewer directions + if (vecSize1 <= vecSize2 ) { min_shell = shell1;} + else { min_shell = shell2;} + if (min_shell.size() > vecSize3){ min_shell = shell3;} + + //most directions + if (vecSize1 >= vecSize2 ) { max_shell = shell1;} + else { max_shell = shell2;} + if (max_shell.size() < vecSize3){ max_shell = shell3;} + + m_MaxDirections = max_shell.size(); + + //SH-order determination + while( ((Interpolation_SHOrder+1)*(Interpolation_SHOrder+2)/2) > min_shell.size() && Interpolation_SHOrder > L ) + Interpolation_SHOrder -= 2 ; + + MITK_INFO << "Interpolation enabeled, using SH of order : " << Interpolation_SHOrder; + + // create target SH-Basis + vnl_matrix * Q = new vnl_matrix(3, max_shell.size()); + ComputeSphericalFromCartesian(Q, max_shell); + + int NumberOfCoeffs = (int)(Interpolation_SHOrder*Interpolation_SHOrder + Interpolation_SHOrder + 2.0)/2.0 + Interpolation_SHOrder; + m_Interpolation_TARGET_SH = new vnl_matrix(max_shell.size(), NumberOfCoeffs); + ComputeSphericalHarmonicsBasis(Q, m_Interpolation_TARGET_SH, Interpolation_SHOrder); + delete Q; + // end creat target SH-Basis + + // create measured-SHBasis + // Shell 1 + vnl_matrix * tempSHBasis; + vnl_matrix_inverse * temp; + + Q = new vnl_matrix(3, shell1.size()); + ComputeSphericalFromCartesian(Q, shell1); + + tempSHBasis = new vnl_matrix(shell1.size(), NumberOfCoeffs); + ComputeSphericalHarmonicsBasis(Q, tempSHBasis, Interpolation_SHOrder); + temp = new vnl_matrix_inverse((*tempSHBasis)); + + m_Interpolation_SHT1_inv = new vnl_matrix(temp->inverse()); + + delete Q; + delete temp; + delete tempSHBasis; + + // Shell 2 + Q = new vnl_matrix(3, shell2.size()); + ComputeSphericalFromCartesian(Q, shell2); + + tempSHBasis = new vnl_matrix(shell2.size(), NumberOfCoeffs); + ComputeSphericalHarmonicsBasis(Q, tempSHBasis, Interpolation_SHOrder); + temp = new vnl_matrix_inverse((*tempSHBasis)); + + m_Interpolation_SHT2_inv = new vnl_matrix(temp->inverse()); + + delete Q; + delete temp; + delete tempSHBasis; + + // Shell 3 + Q = new vnl_matrix(3, shell3.size()); + ComputeSphericalFromCartesian(Q, shell3); + + tempSHBasis = new vnl_matrix(shell3.size(), NumberOfCoeffs); + + ComputeSphericalHarmonicsBasis(Q, tempSHBasis, Interpolation_SHOrder); + + temp = new vnl_matrix_inverse((*tempSHBasis)); + + m_Interpolation_SHT3_inv = new vnl_matrix(temp->inverse()); + + delete Q; + delete temp; + delete tempSHBasis; + + ComputeReconstructionMatrix(max_shell); + return; + }else + { + ComputeReconstructionMatrix(shell1); + } } } if(m_BValueMap.size() > 2 && m_ReconstructionType != Mode_Analytical3Shells) { m_ReconstructionType = Mode_NumericalNShells; } - - this->m_NumberOfBaselineImages = m_BValueMap[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."); + if(m_BValueMap.size() == 2){ + BValueMapIteraotr it = m_BValueMap.begin(); + it++; // skip b0 entry + IndiciesVector shell = it->second; + ComputeReconstructionMatrix(shell); } - this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); } + template< class T, class TG, class TO, int L, int NODF> void DiffusionMultiShellQballReconstructionImageFilter -::BeforeThreadedGenerateData() +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int NumberOfThreads) { itk::TimeProbe clock; clock.Start(); - - if( m_NumberOfGradientDirections < (((L+1)*(L+2))/2) /* && m_NumberOfGradientDirections < 6 */ ) + switch(m_ReconstructionType) { - itkExceptionMacro( << "At least " << ((L+1)*(L+2))/2 << " gradient directions are required" ); + case Mode_Standard1Shell: + StandardOneShellReconstruction(outputRegionForThread); + break; + case Mode_Analytical3Shells: + AnalyticalThreeShellReconstruction(outputRegionForThread); + break; + case Mode_NumericalNShells: + break; } - // 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(); - + MITK_INFO << "Reconstruction in : " << clock.GetTotal() << " s"; } - 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)); + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(ProcessObject::GetOutput(0)); + // Get input gradient image pointer + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( ProcessObject::GetInput(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_BValueMap[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) ); + ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread); + bzeroIterator.GoToBegin(); // Const ImageRegionIterator for input gradient image typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType; GradientIteratorType git(gradientImagePointer, outputRegionForThread ); git.GoToBegin(); + BValueMapIteraotr it = m_BValueMap.begin(); + it++; // skip b0 entry + IndiciesVector SignalIndicies = it->second; + IndiciesVector BZeroIndicies = m_BValueMap[0]; + + int NumbersOfGradientIndicies = SignalIndicies.size(); + 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); + double b0average = 0; + const int b0size = BZeroIndicies.size(); + for(unsigned int i = 0; i SignalVector(m_NumberOfGradientDirections); - if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) + vnl_vector SignalVector(NumbersOfGradientIndicies); + if( (b0average != 0) && (b0average >= 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); + S_S0Normalization(SignalVector, b0average); + Projection1(SignalVector); - // ODF coeffs-vector - vnl_vector coeffs(m_NumberCoefficients); + DoubleLogarithm(SignalVector); // approximate ODF coeffs - coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); + vnl_vector coeffs = ( (*m_CoeffReconstructionMatrix) * SignalVector ); coeffs[0] = 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI)); odf = 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_BValueMap[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) ); - + typedef typename GradientImagesType::PixelType GradientVectorType; + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(ProcessObject::GetOutput(0)); + typename GradientImagesType::Pointer gradientImagePointer = static_cast< GradientImagesType * >( ProcessObject::GetInput(0) ); // Define Image iterators ImageRegionIterator< OutputImageType > odfOutputImageIterator(outputImage, outputRegionForThread); - ImageRegionConstIterator< BZeroImageType > bzeroImageIterator(m_BZeroImage, outputRegionForThread); ImageRegionConstIterator< GradientImagesType > gradientInputImageIterator(gradientImagePointer, outputRegionForThread ); + ImageRegionIterator< BZeroImageType > bzeroIterator(m_BZeroImage, outputRegionForThread); // All iterators seht to Begin of the specific OutputRegion + bzeroIterator.GoToBegin(); odfOutputImageIterator.GoToBegin(); - bzeroImageIterator.GoToBegin(); gradientInputImageIterator.GoToBegin(); // Get Shell Indicies for all non-BZero Gradients // it MUST be a arithmetic progression eg.: 1000, 2000, 3000 BValueMapIteraotr it = m_BValueMap.begin(); it++; // it = b-value = 1000 - IndiciesVector Shell1Indiecies = (*it).second; + IndiciesVector Shell1Indiecies = it->second; it++; // it = b-value = 2000 - IndiciesVector Shell2Indiecies = (*it).second; + IndiciesVector Shell2Indiecies = it->second; it++; // it = b-value = 3000 - IndiciesVector Shell3Indiecies = (*it).second; + IndiciesVector Shell3Indiecies = it->second; + IndiciesVector BZeroIndicies = m_BValueMap[0]; + if(!m_Interpolation_Flag) + { + m_MaxDirections = Shell1Indiecies.size(); + }// else: m_MaxDirection is set in BeforeThreadedGenerateData - // 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 Matrix with E(0) = Shell 1, E(1) = Shell 2, E(2) = Shell 3 + vnl_vector< double > E1(m_MaxDirections); + vnl_vector< double > E2(m_MaxDirections); + vnl_vector< double > E3(m_MaxDirections); - // 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(m_MaxDirections); + vnl_vector BetaValues(m_MaxDirections); + vnl_vector LAValues(m_MaxDirections); + vnl_vector PValues(m_MaxDirections); - 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()); + vnl_vector DataShell1(Shell1Indiecies.size()); + vnl_vector DataShell2(Shell2Indiecies.size()); + vnl_vector DataShell3(Shell3Indiecies.size()); OdfPixelType odf(0.0); + double P2,A,B2,B,P,alpha,beta,lambda, ER1, ER2; + + // iterate overall voxels of the gradient image region while( ! gradientInputImageIterator.IsAtEnd() ) { - if( (bzeroImageIterator.Get() != 0) && (bzeroImageIterator.Get() >= m_Threshold) ) + + GradientVectorType b = gradientInputImageIterator.Get(); + + // calculate for each shell the corresponding b0-averages + double shell1b0Norm =0; + double shell2b0Norm =0; + double shell3b0Norm =0; + double b0average = 0; + const int b0size = BZeroIndicies.size(); + for(unsigned int i = 0; i = b0size / 3 && i < (b0size / 3)*2) shell2b0Norm += b[BZeroIndicies[i]]; + if(i >= (b0size / 3) * 2) shell3b0Norm += b[BZeroIndicies[i]]; + } + shell1b0Norm /= (BZeroIndicies.size()/3); + shell2b0Norm /= (BZeroIndicies.size()/3); + shell3b0Norm /= (BZeroIndicies.size()/3); + b0average = (shell1b0Norm + shell2b0Norm+ shell3b0Norm)/3; + bzeroIterator.Set(b0average); + ++bzeroIterator; + + if( (b0average != 0) && ( b0average >= 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(); + + ///fsl fix --------------------------------------------------- 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]])); - } + DataShell1[i] = static_cast(b[Shell1Indiecies[i]]); + for(int i = 0 ; i < Shell2Indiecies.size(); i++) + DataShell2[i] = static_cast(b[Shell2Indiecies[i]]); + for(int i = 0 ; i < Shell3Indiecies.size(); i++) + DataShell3[i] = static_cast(b[Shell2Indiecies[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(DataShell1, shell1b0Norm); + S_S0Normalization(DataShell2, shell2b0Norm); + S_S0Normalization(DataShell3, shell2b0Norm); + //fsl fix -------------------------------------------ende-- + + /* correct version + for(int i = 0 ; i < Shell1Indiecies.size(); i++) + DataShell1[i] = static_cast(b[Shell1Indiecies[i]]); + for(int i = 0 ; i < Shell2Indiecies.size(); i++) + DataShell2[i] = static_cast(b[Shell2Indiecies[i]]); + for(int i = 0 ; i < Shell3Indiecies.size(); i++) + DataShell3[i] = static_cast(b[Shell3Indiecies[i]]); // Normalize the Signal: Si/S0 - S_S0Normalization(*E,bzeroImageIterator.Get()); + S_S0Normalization(DataShell1, shell1b0Norm); + S_S0Normalization(DataShell2, shell2b0Norm); + S_S0Normalization(DataShell3, shell3b0Norm); + */ - //Implements Eq. [19] and Fig. 4. - Threshold(*E); + if(m_Interpolation_Flag) + { + E1 = ((*m_Interpolation_TARGET_SH) * (*m_Interpolation_SHT1_inv) * (DataShell1)); + E2 = ((*m_Interpolation_TARGET_SH) * (*m_Interpolation_SHT2_inv) * (DataShell2)); + E3 = ((*m_Interpolation_TARGET_SH) * (*m_Interpolation_SHT3_inv) * (DataShell3)); + }else{ + E1 = (DataShell1); + E2 = (DataShell2); + E3 = (DataShell3); + } + //Implements Eq. [19] and Fig. 4. + Projection1(E1); + Projection1(E2); + Projection1(E3); //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; + Projection2(E1,E2,E3); - for( unsigned int i = 0; i< Shell1Indiecies.size(); i++ ) + for( unsigned int i = 0; i< m_MaxDirections; i++ ) { + double e1 = E1.get(i); + double e2 = E2.get(i); + double e3 = E3.get(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; + 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; + PValues.put(i, P); + AlphaValues.put(i, alpha); + BetaValues.put(i, beta); - 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 )); + } + Projection3(PValues, AlphaValues, BetaValues); - PValues->put(i, P); - AlphaValues->put(i, alpha); - BetaValues->put(i, beta); - LAValues->put(i,(lambda * (ER1 < ER2)) + ((1-lambda) * (ER1 >= ER2))); + for(int i = 0 ; i < m_MaxDirections; i++) + { + const double fac = (PValues[i] * 2 ) / (AlphaValues[i] - BetaValues[i]); + lambda = 0.5 + 0.5 * std::sqrt(1 - fac * fac);; + ER1 = std::fabs(lambda * (AlphaValues[i] - BetaValues[i]) + (BetaValues[i] - E1.get(i) )) + + std::fabs(lambda * (AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] - E2.get(i) )) + + std::fabs(lambda * (AlphaValues[i] * AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] * BetaValues[i] - E3.get(i) )); + ER2 = std::fabs((1-lambda) * (AlphaValues[i] - BetaValues[i]) + (BetaValues[i] - E1.get(i) )) + + std::fabs((1-lambda) * (AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] - E2.get(i) )) + + std::fabs((1-lambda) * (AlphaValues[i] * AlphaValues[i] * AlphaValues[i] - BetaValues[i] * BetaValues[i] * BetaValues[i]) + (BetaValues[i] * BetaValues[i] * BetaValues[i] - E3.get(i))); + if(ER1 < ER2) + LAValues.put(i, lambda); + else + LAValues.put(i, 1-lambda); } - Projection2(*PValues, *AlphaValues, *BetaValues); - //Threshold(*AlphaValues); - //Threshold(*BetaValues); + DoubleLogarithm(AlphaValues); + DoubleLogarithm(BetaValues); - DoubleLogarithm(*AlphaValues); - DoubleLogarithm(*BetaValues); + vnl_vector SignalVector(element_product((LAValues) , (AlphaValues)-(BetaValues)) + (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 = element_cast( (*m_ODFSphericalHarmonicBasisMatrix) * coeffs ).data_block(); - odf *= (QBALL_ANAL_RECON_PI*4/NODF); - //Normalize(odf); + odf *= ((QBALL_ANAL_RECON_PI*4)/NODF); } // 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) +void DiffusionMultiShellQballReconstructionImageFilter:: +ComputeSphericalHarmonicsBasis(vnl_matrix * QBallReference, vnl_matrix *SHBasisOutput, int LOrder , vnl_matrix* LaplaciaBaltramiOutput, vnl_vector* SHOrderAssociation, vnl_matrix* SHEigenvalues) { - itk::TimeProbe clock; - GenerateAveragedBZeroImage(outputRegionForThread); + // MITK_INFO << *QBallReference; - 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 k = 0; k <= LOrder; 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); + if(QBallReference){ + 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() +::ComputeReconstructionMatrix(IndiciesVector const & refVector) { 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_BValueMap.begin()); - it++; - const std::vector gradientIndiciesVector = (*it).second; - - int numberOfGradientDirections = gradientIndiciesVector.size(); + int numberOfGradientDirections = refVector.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)); + int NumberOfCoeffs = (int)(LOrder*LOrder + LOrder + 2.0)/2.0 + LOrder; + MITK_INFO << NumberOfCoeffs; + MatrixDoublePtr SHBasisMatrix(new vnl_matrix(numberOfGradientDirections,NumberOfCoeffs)); + SHBasisMatrix->fill(0.0); + VectorIntPtr SHOrderAssociation(new vnl_vector(NumberOfCoeffs)); SHOrderAssociation->fill(0.0); - MatrixDoublePtr LaplacianBaltrami(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + MatrixDoublePtr LaplacianBaltrami(new vnl_matrix(NumberOfCoeffs,NumberOfCoeffs)); LaplacianBaltrami->fill(0.0); - MatrixDoublePtr FRTMatrix(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + MatrixDoublePtr FRTMatrix(new vnl_matrix(NumberOfCoeffs,NumberOfCoeffs)); FRTMatrix->fill(0.0); - MatrixDoublePtr SHEigenvalues(new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients)); + MatrixDoublePtr SHEigenvalues(new vnl_matrix(NumberOfCoeffs,NumberOfCoeffs)); SHEigenvalues->fill(0.0); + MatrixDoublePtr Q(new vnl_matrix(3, numberOfGradientDirections)); + // Convert Cartesian to Spherical Coordinates refVector -> Q + ComputeSphericalFromCartesian(Q.get(), refVector); // SHBasis-Matrix + LaplacianBaltrami-Matrix + SHOrderAssociationVector - ComputeSphericalHarmonicsBasis(Q.get() ,m_SHBasisMatrix, LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); + ComputeSphericalHarmonicsBasis(Q.get() ,SHBasisMatrix.get() , LOrder , LaplacianBaltrami.get(), SHOrderAssociation.get(), SHEigenvalues.get()); // Compute FunkRadon Transformation Matrix Associated to SHBasis Order lj - ComputeFunkRadonTransformationMatrix(SHOrderAssociation.get() ,FRTMatrix.get()); + for(int i=0; i(((m_SHBasisMatrix->transpose()) * (*m_SHBasisMatrix)) + (m_Lambda * (*LaplacianBaltrami)))); + MatrixDoublePtr temp(new vnl_matrix(((SHBasisMatrix->transpose()) * (*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); + MatrixDoublePtr inverse(new vnl_matrix(NumberOfCoeffs,NumberOfCoeffs)); (*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())); + const double factor = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI)); + MatrixDoublePtr SignalReonstructionMatrix (new vnl_matrix((*inverse) * (SHBasisMatrix->transpose()))); + m_CoeffReconstructionMatrix = new vnl_matrix(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*SignalReonstructionMatrix))) )); - m_CoeffReconstructionMatrix = new vnl_matrix(( factor * ((*FRTMatrix) * ((*SHEigenvalues) * (*m_SignalReonstructionMatrix))) )); - - // this code goes to the image adapter coeffs->odfs later + // SH Basis for ODF-reconstruction 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( U->as_matrix() )); + m_ODFSphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,NumberOfCoeffs); + ComputeSphericalHarmonicsBasis(tempPtr.get(), m_ODFSphericalHarmonicBasisMatrix, LOrder); } -template< class T, class TG, class TO, int L, int NODF> -template< class VNLType > -void DiffusionMultiShellQballReconstructionImageFilter -::printMatrix( VNLType * mat ) +template< class T, class TG, class TO, int L, int NOdfDirections> +void DiffusionMultiShellQballReconstructionImageFilter +::ComputeSphericalFromCartesian(vnl_matrix * Q, IndiciesVector const & refShell) { - - std::stringstream stream; - - for(int i = 0 ; i < mat->rows(); i++) + for(int i = 0; i < refShell.size(); i++) { - stream.str(""); - for(int j = 0; j < mat->cols(); j++) - { - stream << (*mat)(i,j) << " "; - } - + double x = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(0); + double y = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(1); + double z = m_GradientDirectionContainer->ElementAt(refShell[i]).normalize().get(2); + double cart[3]; + mitk::sh::Cart2Sph(x,y,z,cart); + (*Q)(0,i) = cart[0]; + (*Q)(1,i) = cart[1]; + (*Q)(2,i) = cart[2]; } - MITK_INFO << stream.str(); } + template< class T, class TG, class TO, int L, int NODF> bool DiffusionMultiShellQballReconstructionImageFilter ::CheckDuplicateDiffusionGradients() { bool value = false; BValueMapIteraotr mapIterator = m_BValueMap.begin(); + mapIterator++; while(mapIterator != m_BValueMap.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; } +// corresponding directions between shells (e.g. dir1_shell1 vs dir1_shell2) differ more than 1 degree. +template< class T, class TG, class TO, int L, int NODF> +bool DiffusionMultiShellQballReconstructionImageFilter +::CheckForDifferingShellDirections() +{ + bool interp_flag = false; + + BValueMapIteraotr mapIterator = m_BValueMap.begin(); + mapIterator++; + std::vector shell1 = mapIterator->second; + mapIterator++; + std::vector shell2 = mapIterator->second; + mapIterator++; + std::vector shell3 = mapIterator->second; + + for (int i=0; i< shell1.size(); i++) + if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell1[i]), m_GradientDirectionContainer->ElementAt(shell2[i]))) <= 0.9998) {interp_flag=true; break;} + for (int i=0; i< shell1.size(); i++) + if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell1[i]), m_GradientDirectionContainer->ElementAt(shell3[i]))) <= 0.9998) {interp_flag=true; break;} + for (int i=0; i< shell1.size(); i++) + if (fabs(dot(m_GradientDirectionContainer->ElementAt(shell2[i]), m_GradientDirectionContainer->ElementAt(shell3[i]))) <= 0.9998) {interp_flag=true; break;} + return interp_flag; +} + + 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 index fb7a499f38..a54d102e7e 100644 --- a/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,229 +1,240 @@ /*=================================================================== 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" #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 > BValueMap; typedef std::map >::iterator BValueMapIteraotr; 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 S_S0Normalization( vnl_vector & vec, double 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); + void Projection1(vnl_vector & vec, 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); + void Projection2( vnl_vector & E1, vnl_vector & E2, vnl_vector & E3, double delta = 0.01); + void Projection3( 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); void SetBValueMap(BValueMap map){this->m_BValueMap = map;} protected: DiffusionMultiShellQballReconstructionImageFilter(); - ~DiffusionMultiShellQballReconstructionImageFilter() {}; + ~DiffusionMultiShellQballReconstructionImageFilter() { }; void PrintSelf(std::ostream& os, Indent indent) const; - void ComputeReconstructionMatrix(); + void ComputeReconstructionMatrix(IndiciesVector const & refVector); + void ComputeODFSHBasis(); 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(); + bool CheckForDifferingShellDirections(); + void ComputeSphericalHarmonicsBasis(vnl_matrix* QBallReference, vnl_matrix* SHBasisOutput, int Lorder , vnl_matrix* LaplaciaBaltramiOutput =0 , vnl_vector* SHOrderAssociation =0 , vnl_matrix * SHEigenvalues =0); + //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 }; + // Interpolation + bool m_Interpolation_Flag; + CoefficientMatrixType m_Interpolation_SHT1_inv; + CoefficientMatrixType m_Interpolation_SHT2_inv; + CoefficientMatrixType m_Interpolation_SHT3_inv; + CoefficientMatrixType m_Interpolation_TARGET_SH; + int m_MaxDirections; //CoefficientMatrixType m_ReconstructionMatrix; CoefficientMatrixType m_CoeffReconstructionMatrix; CoefficientMatrixType m_ODFSphericalHarmonicBasisMatrix; - CoefficientMatrixType m_SignalReonstructionMatrix; - CoefficientMatrixType m_SHBasisMatrix; + //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; BValueMap m_BValueMap; double m_Lambda; bool m_IsHemisphericalArrangementOfGradientDirections; bool m_IsArithmeticProgession; - int m_NumberCoefficients; + //int m_NumberCoefficients; ReconstructionType m_ReconstructionType; - - template< class VNLType > - void printMatrix( VNLType * mat ); - //------------------------- VNL-function ------------------------------------ - template - static vnl_vector< WntValue> element_cast (vnl_vector< CurrentValue> const& v1) + vnl_vector< WntValue> element_cast (vnl_vector< CurrentValue> const& v1) { vnl_vector result(v1.size()); for(int i = 0 ; i < v1.size(); i++) result[i] = static_cast< WntValue>(v1[i]); return result; } + template + 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 ; + } + + void ComputeSphericalFromCartesian(vnl_matrix * Q, const IndiciesVector & refShell); + + }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_