itkDiffusionQballPrepareVisualizationImageFilter.cpp
Go to the documentation of this file.
00001 00002 #ifndef __itkDiffusionQballPrepareVisualizationImageFilter_cpp 00003 #define __itkDiffusionQballPrepareVisualizationImageFilter_cpp 00004 00005 #include <time.h> 00006 #include <stdio.h> 00007 #include <stdlib.h> 00008 00009 #include "itkDiffusionQballPrepareVisualizationImageFilter.h" 00010 #include "itkImageRegionConstIterator.h" 00011 #include "itkImageRegionConstIteratorWithIndex.h" 00012 #include "itkImageRegionIterator.h" 00013 #include "itkArray.h" 00014 #include "vnl/vnl_vector.h" 00015 #include "itkOrientationDistributionFunction.h" 00016 00017 namespace itk { 00018 00019 //#define QBALL_RECON_PI 3.14159265358979323846 00020 00021 template< class TOdfPixelType, int NrOdfDirections> 00022 DiffusionQballPrepareVisualizationImageFilter< TOdfPixelType, NrOdfDirections> 00023 ::DiffusionQballPrepareVisualizationImageFilter() : 00024 m_Threshold(0), 00025 m_ScaleByGfaType(GfaFilterType::GFA_STANDARD), 00026 m_DoScaleGfa(false), 00027 m_GfaParam1(2), 00028 m_GfaParam2(1) 00029 { 00030 // At least 1 inputs is necessary for a vector image. 00031 // For images added one at a time we need at least six 00032 this->SetNumberOfRequiredInputs( 1 ); 00033 } 00034 00035 template< class TOdfPixelType, 00036 int NrOdfDirections> 00037 void DiffusionQballPrepareVisualizationImageFilter< TOdfPixelType, 00038 NrOdfDirections> 00039 ::BeforeThreadedGenerateData() 00040 { 00041 if( m_NormalizationMethod == PV_GLOBAL_MAX ) 00042 { 00043 typename InputImageType::Pointer inputImagePointer = NULL; 00044 inputImagePointer = static_cast< InputImageType * >( 00045 this->ProcessObject::GetInput(0) ); 00046 00047 typename GfaFilterType::Pointer filter = GfaFilterType::New(); 00048 filter->SetInput(inputImagePointer); 00049 filter->SetNumberOfThreads(4); 00050 filter->SetComputationMethod(GfaFilterType::GFA_MAX_ODF_VALUE); 00051 filter->Update(); 00052 00053 typedef typename itk::MinimumMaximumImageCalculator< typename GfaFilterType::OutputImageType > 00054 MaxFilterType; 00055 typename MaxFilterType::Pointer maxFilter = MaxFilterType::New(); 00056 maxFilter->SetImage(filter->GetOutput()); 00057 maxFilter->ComputeMaximum(); 00058 00059 m_GlobalInputMaximum = maxFilter->GetMaximum(); 00060 } 00061 00062 //if(m_DoScaleGfa) 00063 { 00064 typename InputImageType::Pointer inputImagePointer = NULL; 00065 inputImagePointer = static_cast< InputImageType * >( 00066 this->ProcessObject::GetInput(0) ); 00067 00068 typename GfaFilterType::Pointer filter = GfaFilterType::New(); 00069 filter->SetInput(inputImagePointer); 00070 filter->SetNumberOfThreads(4); 00071 filter->SetComputationMethod(m_ScaleByGfaType); 00072 filter->SetParam1(m_GfaParam1); 00073 filter->SetParam2(m_GfaParam2); 00074 filter->Update(); 00075 m_GfaImage = filter->GetOutput(); 00076 } 00077 } 00078 00079 template< class TOdfPixelType, 00080 int NrOdfDirections> 00081 void DiffusionQballPrepareVisualizationImageFilter< TOdfPixelType, 00082 NrOdfDirections> 00083 ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, 00084 int ) 00085 { 00086 typename OutputImageType::Pointer outputImage = 00087 static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); 00088 ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); 00089 oit.GoToBegin(); 00090 00091 typedef itk::OrientationDistributionFunction<TOdfPixelType,NrOdfDirections> OdfType; 00092 typedef ImageRegionConstIterator< InputImageType > InputIteratorType; 00093 typedef typename InputImageType::PixelType OdfVectorType; 00094 typename InputImageType::Pointer inputImagePointer = NULL; 00095 inputImagePointer = static_cast< InputImageType * >( 00096 this->ProcessObject::GetInput(0) ); 00097 InputIteratorType git(inputImagePointer, outputRegionForThread ); 00098 git.GoToBegin(); 00099 00100 typedef ImageRegionConstIterator< GfaImageType > GfaIteratorType; 00101 GfaIteratorType gfaIt(m_GfaImage, outputRegionForThread); 00102 00103 while( !git.IsAtEnd() ) 00104 { 00105 OdfVectorType b = git.Get(); 00106 OdfType odf = b.GetDataPointer(); 00107 00108 switch( m_NormalizationMethod ) 00109 { 00110 case PV_NONE: 00111 { 00112 break; 00113 } 00114 case PV_MAX: 00115 { 00116 odf = odf.MaxNormalize(); 00117 break; 00118 } 00119 case PV_MIN_MAX: 00120 { 00121 odf = odf.MinMaxNormalize(); 00122 break; 00123 } 00124 case PV_GLOBAL_MAX: 00125 { 00126 odf *= 1.0/m_GlobalInputMaximum; 00127 break; 00128 } 00129 case PV_MIN_MAX_INVERT: 00130 { 00131 odf = odf.MinMaxNormalize(); 00132 for(int i=0; i<NrOdfDirections; i++) 00133 { 00134 odf[i] = 1.0 - odf[i]; 00135 } 00136 break; 00137 } 00138 } 00139 00140 if(m_DoScaleGfa) 00141 { 00142 odf *= gfaIt.Get(); 00143 ++gfaIt; 00144 } 00145 00146 odf *= 0.5; 00147 oit.Set( odf.GetDataPointer() ); 00148 ++oit; 00149 ++git; // Gradient image iterator 00150 } 00151 00152 std::cout << "One Thread finished extraction" << std::endl; 00153 } 00154 00155 template< class TOdfPixelType, 00156 int NrOdfDirections> 00157 void DiffusionQballPrepareVisualizationImageFilter< TOdfPixelType, 00158 NrOdfDirections> 00159 ::PrintSelf(std::ostream& os, Indent indent) const 00160 { 00161 Superclass::PrintSelf(os,indent); 00162 os << indent << "m_Threshold: " << 00163 m_Threshold << std::endl; 00164 } 00165 00166 } 00167 00168 #endif // __itkDiffusionQballPrepareVisualizationImageFilter_cpp
