1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
16 #ifndef __itkDiffusionQballReconstructionImageFilter_txx
17 #define __itkDiffusionQballReconstructionImageFilter_txx
19 #include "itkDiffusionQballReconstructionImageFilter.h"
20 #include "itkImageRegionConstIterator.h"
21 #include "itkImageRegionConstIteratorWithIndex.h"
22 #include "itkImageRegionIterator.h"
24 #include "vnl/vnl_vector.h"
25 #include "itkPointShell.h"
27 #define _USE_MATH_DEFINES
32 #define QBALL_RECON_PI M_PI
34 template< class TReferenceImagePixelType,
35 class TGradientImagePixelType,
38 int NrBasisFunctionCenters>
39 DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
40 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
41 NrBasisFunctionCenters>
42 ::DiffusionQballReconstructionImageFilter() :
43 m_GradientDirectionContainer(NULL),
44 m_NumberOfGradientDirections(0),
45 m_NumberOfEquatorSamplingPoints(0),
46 m_NumberOfBaselineImages(1),
47 m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
49 m_GradientImageTypeEnumeration(Else),
50 m_DirectionsDuplicated(false)
52 // At least 1 inputs is necessary for a vector image.
53 // For images added one at a time we need at least six
54 this->SetNumberOfRequiredInputs( 1 );
57 template< class TReferenceImagePixelType,
58 class TGradientImagePixelType,
61 int NrBasisFunctionCenters>
62 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
63 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
64 NrBasisFunctionCenters>
65 ::BeforeThreadedGenerateData()
67 // If we have more than 2 inputs, then each input, except the first is a
68 // gradient image. The number of gradient images must match the number of
69 // gradient directions.
70 const unsigned int numberOfInputs = this->GetNumberOfInputs();
72 // There need to be at least 6 gradient directions to be able to compute the
74 if( m_NumberOfGradientDirections < 1 )
76 itkExceptionMacro( << "Your image contains no diffusion gradients!" );
79 // If there is only 1 gradient image, it must be an itk::VectorImage. Otherwise
80 // we must have a container of (numberOfInputs-1) itk::Image. Check to make sure
81 if ( numberOfInputs == 1
82 && m_GradientImageTypeEnumeration != GradientIsInASingleImage )
84 std::string gradientImageClassName(
85 this->ProcessObject::GetInput(0)->GetNameOfClass());
86 if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
89 "There is only one Gradient image. I expect that to be a VectorImage. "
90 << "But its of type: " << gradientImageClassName );
94 // Compute reconstruction matrix that is multiplied to the data-vector
95 // each voxel in order to reconstruct the ODFs
96 this->ComputeReconstructionMatrix();
98 // Allocate the b-zero image
99 m_BZeroImage = BZeroImageType::New();
100 if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
102 typename ReferenceImageType::Pointer img = static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0));
103 m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
104 m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin
105 m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction
106 m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
107 m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
109 // The gradients are specified in a single multi-component image
110 else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage )
112 typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
113 this->ProcessObject::GetInput(0) );
114 m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
115 m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin
116 m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction
117 m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
118 m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
121 m_BZeroImage->Allocate();
125 template< class TReferenceImagePixelType,
126 class TGradientImagePixelType,
129 int NrBasisFunctionCenters>
130 typename itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::OdfPixelType
131 itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::Normalize( OdfPixelType odf, typename NumericTraits<ReferencePixelType>::AccumulateType b0 )
133 switch( m_NormalizationMethod )
136 // divide by sum to retreive a PDF
144 case QBR_B_ZERO_B_VALUE:
146 for(int i=0; i<NrOdfDirections; i++)
148 odf[i] = ((TOdfPixelType)log((TOdfPixelType)b0)-odf[i])/m_BValue;
153 // divide by b-zero value of voxel
160 // no normalization of ODF
172 template< class TReferenceImagePixelType,
173 class TGradientImagePixelType,
176 int NrBasisFunctionCenters >
177 vnl_vector<TOdfPixelType> itk::DiffusionQballReconstructionImageFilter<TReferenceImagePixelType, TGradientImagePixelType, TOdfPixelType, NrOdfDirections, NrBasisFunctionCenters>::PreNormalize( vnl_vector<TOdfPixelType> vec )
179 switch( m_NormalizationMethod )
181 // standard: no normalization before reconstruction
188 case QBR_B_ZERO_B_VALUE:
191 for(int i=0; i<n; i++)
193 vec[i] = log(vec[i]);
198 // no normalization before reconstruction here
204 // no normalization before reconstruction here
215 template< class TReferenceImagePixelType,
216 class TGradientImagePixelType,
219 int NrBasisFunctionCenters>
220 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
221 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
222 NrBasisFunctionCenters>
223 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
226 // init output and b-zero iterators
227 typename OutputImageType::Pointer outputImage =
228 static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput());
229 ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
231 ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
233 vnl_vector<TOdfPixelType> B(m_NumberOfGradientDirections);
236 // 1. Gradients specified in multiple images
237 // 'n' iterators for each of the gradient images
238 // 2. Gradients specified in a single multi-component image
239 // one iterator for all gradient directions
241 if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
243 // b-zero reference image iterator
244 ImageRegionConstIterator< ReferenceImageType >
245 it(static_cast< ReferenceImageType * >(this->ProcessObject::GetInput(0)),
246 outputRegionForThread);
249 // fill vector with gradient iterators
250 typedef ImageRegionConstIterator< GradientImageType > GradientIteratorType;
251 std::vector< GradientIteratorType * > gradientItContainer;
252 for( unsigned int i = 1; i<= m_NumberOfGradientDirections; i++ )
254 // adapt index in case we have a duplicated shell
256 if(m_DirectionsDuplicated)
257 index = i % (m_NumberOfGradientDirections/2);
259 // init and pushback current input image iterator
260 typename GradientImageType::Pointer gradientImagePointer = NULL;
261 // dynamic_cast would be nice, static because of SGI
262 gradientImagePointer = static_cast< GradientImageType * >(
263 this->ProcessObject::GetInput(index) );
264 GradientIteratorType *git = new GradientIteratorType(
265 gradientImagePointer, outputRegionForThread );
267 gradientItContainer.push_back(git);
270 // Following loop does the actual reconstruction work in each voxel
271 // (Tuch, Q-Ball Reconstruction [1])
272 while( !it.IsAtEnd() )
275 // b-zero reference value
276 ReferencePixelType b0 = it.Get();
279 OdfPixelType odf(0.0);
281 // threshold on reference value to suppress noisy regions
282 if( (b0 != 0) && (b0 >= m_Threshold) )
285 // fill array of diffusion measurements
286 for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
288 GradientPixelType b = gradientItContainer[i]->Get();
289 B[i] = static_cast<TOdfPixelType>(b);
290 ++(*gradientItContainer[i]);
293 // pre-normalization according to m_NormalizationMethod
296 // actual reconstruction
297 odf = ( (*m_ReconstructionMatrix) * B ).data_block();
299 // post-normalization according to m_NormalizationMethod
304 // in case we fall below threshold, we just increment to next voxel
305 for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
307 ++(*gradientItContainer[i]);
311 for (unsigned int i=0; i<odf.Size(); i++)
312 if (odf.GetElement(i)!=odf.GetElement(i))
315 // set and increment output iterators
324 for( unsigned int i = 0; i< gradientItContainer.size(); i++ )
326 delete gradientItContainer[i];
329 // The gradients are specified in a single multi-component image
330 else if( m_GradientImageTypeEnumeration == GradientIsInASingleImage )
333 // init input iterator
334 typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
335 typedef typename GradientImagesType::PixelType GradientVectorType;
336 typename GradientImagesType::Pointer gradientImagePointer = NULL;
337 // dynamic_cast would be nice, static because of SGI
338 gradientImagePointer = static_cast< GradientImagesType * >(
339 this->ProcessObject::GetInput(0) );
340 GradientIteratorType git(gradientImagePointer, outputRegionForThread );
343 // set of indicies each for the baseline images and gradient images
344 std::vector<unsigned int> baselineind; // contains baseline indicies
345 std::vector<unsigned int> gradientind; // contains gradient indicies
346 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
347 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
349 if(gdcit.Value().one_norm() <= 0.0)
351 baselineind.push_back(gdcit.Index());
355 gradientind.push_back(gdcit.Index());
359 // in case we have a duplicated shell, we also duplicate or indices
360 if( m_DirectionsDuplicated )
362 int gradIndSize = gradientind.size();
363 for(int i=0; i<gradIndSize; i++)
364 gradientind.push_back(gradientind[i]);
367 // Following loop does the actual reconstruction work in each voxel
368 // (Tuch, Q-Ball Reconstruction [1])
369 while( !git.IsAtEnd() )
371 // current vector of diffusion measurements
372 GradientVectorType b = git.Get();
374 // average of current b-zero reference values
375 typename NumericTraits<ReferencePixelType>::AccumulateType b0 = NumericTraits<ReferencePixelType>::Zero;
376 for(unsigned int i = 0; i < baselineind.size(); ++i)
378 b0 += b[baselineind[i]];
380 b0 /= this->m_NumberOfBaselineImages;
382 // init resulting ODF
383 OdfPixelType odf(0.0);
385 // threshold on reference value to suppress noisy regions
386 if( (b0 != 0) && (b0 >= m_Threshold) )
388 for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
390 B[i] = static_cast<TOdfPixelType>(b[gradientind[i]]);
393 // pre-normalization according to m_NormalizationMethod
396 // actual reconstruction
397 odf = ( (*m_ReconstructionMatrix) * B ).data_block();
399 // post-normalization according to m_NormalizationMethod
400 odf = Normalize(odf, b0);
404 for (unsigned int i=0; i<odf.Size(); i++)
405 if (odf.GetElement(i)!=odf.GetElement(i))
408 // set and increment output iterators
413 ++git; // Gradient image iterator
417 std::cout << "One Thread finished reconstruction" << std::endl;
420 template< class TReferenceImagePixelType,
421 class TGradientImagePixelType,
424 int NrBasisFunctionCenters>
425 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
426 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
427 NrBasisFunctionCenters>
428 ::ComputeReconstructionMatrix()
431 if( m_NumberOfGradientDirections < 1 )
433 itkExceptionMacro( << "Your image contains no diffusion gradients!" );
437 // check for duplicate diffusion gradients
438 bool warning = false;
439 for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
440 gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
442 for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
443 gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
445 if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
447 itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
455 // handle acquisition schemes where only half of the spherical
456 // shell is sampled by the gradient directions. In this case,
457 // each gradient direction is duplicated in negative direction.
458 vnl_vector<double> centerMass(3);
459 centerMass.fill(0.0);
461 for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
462 gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
464 if(gdcit1.Value().one_norm() > 0.0)
466 centerMass += gdcit1.Value();
471 if(centerMass.two_norm() > 0.1)
473 m_DirectionsDuplicated = true;
474 m_NumberOfGradientDirections *= 2;
478 // set default number of equator sampling points if needed
479 if(!this->m_NumberOfEquatorSamplingPoints)
480 this->m_NumberOfEquatorSamplingPoints
481 = (int) ceil((double)sqrt(8*QBALL_RECON_PI*this->m_NumberOfGradientDirections));
483 vnl_matrix<double>* Q =
484 new vnl_matrix<double>(3, m_NumberOfGradientDirections);
487 // Fill matrix Q with gradient directions
489 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
490 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
492 if(gdcit.Value().one_norm() > 0.0)
494 (*Q)(0,i) = gdcit.Value().get(0);
495 (*Q)(1,i) = gdcit.Value().get(1);
496 (*Q)(2,i++) = gdcit.Value().get(2);
499 if(m_DirectionsDuplicated)
501 for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
502 gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
504 if(gdcit.Value().one_norm() > 0.0)
506 (*Q)(0,i) = -gdcit.Value().get(0);
507 (*Q)(1,i) = -gdcit.Value().get(1);
508 (*Q)(2,i++) = -gdcit.Value().get(2);
515 vnl_matrix_fixed<double, 3, NOdfDirections>* U =
516 itk::PointShell<NOdfDirections, vnl_matrix_fixed<double, 3, NOdfDirections> >::DistributePointShell();
518 vnl_matrix_fixed<double, 3, NBasisFunctionCenters>* V =
519 itk::PointShell<NBasisFunctionCenters, vnl_matrix_fixed<double, 3, NBasisFunctionCenters> >::DistributePointShell();
521 // calculate sampling points on the equator perpendicular to z-axis
522 vnl_matrix<double> *C
523 = new vnl_matrix<double>(3, m_NumberOfEquatorSamplingPoints);
524 for(unsigned int i=0; i<m_NumberOfEquatorSamplingPoints; i++)
526 double theta = i * (2*QBALL_RECON_PI / m_NumberOfEquatorSamplingPoints);
527 (*C)(0,i) = cos(theta);
528 (*C)(1,i) = sin(theta);
529 (*C)(2,i) = NumericTraits<double>::Zero;
532 // rotate the sampling points to each directions of interest
533 vnl_matrix<double> *S
534 = new vnl_matrix<double>(3,m_NumberOfEquatorSamplingPoints*NOdfDirections);
537 vnl_vector_fixed<double,3> z(NumericTraits<double>::Zero);
538 z.put(2,NumericTraits<double>::One);
539 vnl_matrix_fixed<double,3,3> eye(NumericTraits<double>::Zero);
540 eye.fill_diagonal(NumericTraits<double>::One);
541 for(int i=0; i<NOdfDirections; i++)
543 vnl_vector_fixed<double,3> ui = (*U).get_column(i);
544 vnl_matrix<double> *RC
545 = new vnl_matrix<double>(3,m_NumberOfEquatorSamplingPoints);
547 if( (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1) != 0 )
549 vnl_matrix_fixed<double,3,3> R;
550 R.set_column(0, (z+ui)*(z(0)+ui(0)));
551 R.set_column(1, (z+ui)*(z(1)+ui(1)));
552 R.set_column(2, (z+ui)*(z(2)+ui(2)));
553 R /= (z(0)*ui(0)+z(1)*ui(1)+z(2)*ui(2)+1);
561 (*S).set_columns(i*m_NumberOfEquatorSamplingPoints, *RC);
565 // determine interpolation kernel width first
566 // use to determine diffusion measurement contribution to each of the kernels
567 vnl_matrix<double> *H_plus
568 = new vnl_matrix<double>(NBasisFunctionCenters,m_NumberOfGradientDirections);
570 double maxSigma = QBALL_RECON_PI/6;
571 double bestSigma = maxSigma;
574 double stepsize = 0.01;
576 double minCondition = NumericTraits<double>::max();
578 vnl_matrix<double> *H
579 = new vnl_matrix<double>(m_NumberOfGradientDirections,NBasisFunctionCenters,(double)0);
584 for( double sigma=start; sigma<maxSigma; sigma+=stepsize)
586 vnl_matrix<double> *tmpH
587 = new vnl_matrix<double>(m_NumberOfGradientDirections,NBasisFunctionCenters);
589 for(unsigned int r=0; r<m_NumberOfGradientDirections; r++)
591 for(int c=0; c<NBasisFunctionCenters; c++)
593 double qtv = (*Q)(0,r)*(*V)(0,c)
594 + (*Q)(1,r)*(*V)(1,c)
595 + (*Q)(2,r)*(*V)(2,c);
596 qtv = (qtv<-1.0) ? -1.0 : ( (qtv>1.0) ? 1.0 : qtv);
597 double x = acos(qtv);
598 (*tmpH)(r,c) = (1.0/(sigma*sqrt(2.0*QBALL_RECON_PI)))
599 *exp((-x*x)/(2*sigma*sigma));
603 vnl_svd<double> *solver;
604 if(m_NumberOfGradientDirections>NBasisFunctionCenters)
606 solver = new vnl_svd<double>(*tmpH);
610 solver = new vnl_svd<double>(tmpH->transpose());
612 double condition = solver->sigma_max() / solver->sigma_min();
614 std::cout << sigma << ": " << condition << std::endl;
616 if( condition < minCondition )
618 minCondition = condition;
624 // optimum assumed to be hit after condition increased 3 times
625 if (++increasing>3) break;
631 vnl_matrix_inverse<double> *pseudoInverse
632 = new vnl_matrix_inverse<double>(*H);
633 (*H_plus) = pseudoInverse->pinverse();
635 std::cout << "choosing sigma = " << bestSigma << std::endl;
639 // this is the contribution of each kernel to each sampling point on the
641 vnl_matrix<double> *G
642 = new vnl_matrix<double>(m_NumberOfEquatorSamplingPoints*NOdfDirections,NBasisFunctionCenters);
645 for(unsigned int r=0; r<m_NumberOfEquatorSamplingPoints*NOdfDirections; r++)
647 for(int c=0; c<NBasisFunctionCenters; c++)
649 double stv = (*S)(0,r)*(*V)(0,c)
650 + (*S)(1,r)*(*V)(1,c)
651 + (*S)(2,r)*(*V)(2,c);
652 stv = (stv<-1.0) ? -1.0 : ( (stv>1.0) ? 1.0 : stv);
653 double x = acos(stv);
654 (*G)(r,c) = (1.0/(bestSigma*sqrt(2.0*QBALL_RECON_PI)))
655 *exp((-x*x)/(2*bestSigma*bestSigma));
660 vnl_matrix<double> *GH_plus =
661 new vnl_matrix<double>(m_NumberOfEquatorSamplingPoints*NOdfDirections,m_NumberOfGradientDirections);
663 // simple matrix multiplication, manual cause of stack overflow using operator
664 for (unsigned i = 0; i < m_NumberOfEquatorSamplingPoints*NOdfDirections; ++i)
665 for (unsigned j = 0; j < m_NumberOfGradientDirections; ++j)
667 double accum = (*G)(i,0) * (*H_plus)(0,j);
668 for (unsigned k = 1; k < NOdfDirections; ++k)
669 accum += (*G)(i,k) * (*H_plus)(k,j);
670 (*GH_plus)(i,j) = accum;
673 typename vnl_matrix<double>::iterator it3;
674 for( it3 = (*GH_plus).begin(); it3 != (*GH_plus).end(); it3++)
680 // this is an addition to the original tuch algorithm
681 for(unsigned int i=0; i<NOdfDirections*m_NumberOfEquatorSamplingPoints; i++)
683 vnl_vector< double > r = GH_plus->get_row(i);
685 GH_plus->set_row(i,r);
688 m_ReconstructionMatrix
689 = new vnl_matrix<TOdfPixelType>(NOdfDirections,m_NumberOfGradientDirections,0.0);
690 for(int i=0; i<NOdfDirections; i++)
692 for(unsigned int j=0; j<m_NumberOfGradientDirections; j++)
694 for(unsigned int k=0; k<m_NumberOfEquatorSamplingPoints; k++)
696 (*m_ReconstructionMatrix)(i,j) += (TOdfPixelType)(*GH_plus)(m_NumberOfEquatorSamplingPoints*i+k,j);
701 // this is also an addition to the original tuch algorithm
702 for(int i=0; i<NOdfDirections; i++)
704 vnl_vector< TOdfPixelType > r = m_ReconstructionMatrix->get_row(i);
706 m_ReconstructionMatrix->set_row(i,r);
708 std::cout << "Reconstruction Matrix computed." << std::endl;
712 template< class TReferenceImagePixelType,
713 class TGradientImagePixelType,
716 int NrBasisFunctionCenters>
717 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
718 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
719 NrBasisFunctionCenters>
720 ::AddGradientImage( const GradientDirectionType &gradientDirection,
721 const GradientImageType *gradientImage )
723 // Make sure crazy users did not call both AddGradientImage and
725 if( m_GradientImageTypeEnumeration == GradientIsInASingleImage)
727 itkExceptionMacro( << "Cannot call both methods:"
728 << "AddGradientImage and SetGradientImage. Please call only one of them.");
731 // If the container to hold the gradient directions hasn't been allocated
733 if( !this->m_GradientDirectionContainer )
735 this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
738 this->m_NumberOfGradientDirections = m_GradientDirectionContainer->Size();
740 m_GradientDirectionContainer->InsertElement( this->m_NumberOfGradientDirections,
741 gradientDirection / gradientDirection.two_norm() );
743 this->ProcessObject::SetNthInput( this->m_NumberOfGradientDirections,
744 const_cast< GradientImageType* >(gradientImage) );
746 m_GradientImageTypeEnumeration = GradientIsInManyImages;
749 template< class TReferenceImagePixelType,
750 class TGradientImagePixelType,
753 int NrBasisFunctionCenters>
754 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
755 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
756 NrBasisFunctionCenters>
757 ::SetGradientImage(const GradientDirectionContainerType *gradientDirection,
758 const GradientImagesType *gradientImage )
760 // Make sure crazy users did not call both AddGradientImage and
762 if( m_GradientImageTypeEnumeration == GradientIsInManyImages )
764 itkExceptionMacro( << "Cannot call both methods:"
765 << "AddGradientImage and SetGradientImage. Please call only one of them.");
768 this->m_GradientDirectionContainer = GradientDirectionContainerType::New();
769 for(GradientDirectionContainerType::ConstIterator it = gradientDirection->Begin();
770 it != gradientDirection->End(); it++)
772 this->m_GradientDirectionContainer->push_back(it.Value());
775 unsigned int numImages = gradientDirection->Size();
776 this->m_NumberOfBaselineImages = 0;
777 for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
778 it != this->m_GradientDirectionContainer->End(); it++)
780 if(it.Value().one_norm() <= 0.0)
782 this->m_NumberOfBaselineImages++;
784 else // Normalize non-zero gradient directions
786 it.Value() = it.Value() / it.Value().two_norm();
790 this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
792 // ensure that the gradient image we received has as many components as
793 // the number of gradient directions
794 if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
796 itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
797 << "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
798 << " directions specified but image has " << gradientImage->GetVectorLength()
802 this->ProcessObject::SetNthInput( 0,
803 const_cast< GradientImagesType* >(gradientImage) );
804 m_GradientImageTypeEnumeration = GradientIsInASingleImage;
808 template< class TReferenceImagePixelType,
809 class TGradientImagePixelType,
812 int NrBasisFunctionCenters>
813 void DiffusionQballReconstructionImageFilter< TReferenceImagePixelType,
814 TGradientImagePixelType, TOdfPixelType, NrOdfDirections,
815 NrBasisFunctionCenters>
816 ::PrintSelf(std::ostream& os, Indent indent) const
819 std::locale originalLocale = os.getloc();
822 Superclass::PrintSelf(os,indent);
824 os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
825 if ( m_GradientDirectionContainer )
827 os << indent << "GradientDirectionContainer: "
828 << m_GradientDirectionContainer << std::endl;
833 "GradientDirectionContainer: (Gradient directions not set)" << std::endl;
835 os << indent << "NumberOfGradientDirections: " <<
836 m_NumberOfGradientDirections << std::endl;
837 os << indent << "NumberOfBaselineImages: " <<
838 m_NumberOfBaselineImages << std::endl;
839 os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
840 os << indent << "BValue: " << m_BValue << std::endl;
841 if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages )
843 os << indent << "Gradient images haven been supplied " << std::endl;
845 else if ( this->m_GradientImageTypeEnumeration == GradientIsInManyImages )
847 os << indent << "A multicomponent gradient image has been supplied" << std::endl;
850 os.imbue( originalLocale );
855 #endif // __itkDiffusionQballReconstructionImageFilter_txx