DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
itkProfileFromSPFImageFilter.hxx
Go to the documentation of this file.
1 
19 #ifndef __itkProfileFromSPFImageFilter_hxx
20 #define __itkProfileFromSPFImageFilter_hxx
21 
22 #include <itkProgressReporter.h>
25 #include "utl.h"
27 
28 namespace itk
29 {
30 
31 template< class TInputImage, class TOutputImage >
32 void
35 {
36  utlShowPosition(this->GetDebug());
37  Superclass::GenerateOutputInformation();
38  typename TOutputImage::Pointer outputPtr = this->GetOutput();
39  unsigned int numberOfComponentsPerPixel = 0;
40  if (this->m_Orientations->Rows()==0)
41  {
42  // output SH coefficients
43  numberOfComponentsPerPixel = utl::RankToDimSH(this->GetSHRank());
44  }
45  else
46  {
47  // output samples
48  numberOfComponentsPerPixel = this->m_Orientations->Rows();
49  utlGlobalException(m_RadiusVector->size()>0 && m_RadiusVector->size()!=this->m_Orientations->Rows(), "wrong size of m_RadiusVector or m_Orientations");
50  }
51  outputPtr->SetNumberOfComponentsPerPixel(numberOfComponentsPerPixel);
52 }
53 
54 template< class TInputImage, class TOutputImage >
55 void
58 {
59  utlShowPosition(this->GetDebug());
60  Superclass::VerifyInputParameters();
61  utlGlobalException(this->m_Radius<0 && this->m_RadiusVector->size()==0, "need to set radius or radiusVector");
62 }
63 
64 template< class TInputImage, class TOutputImage >
65 typename LightObject::Pointer
68 {
69  typename LightObject::Pointer loPtr = Superclass::InternalClone();
70 
71  typename Self::Pointer rval = dynamic_cast<Self *>(loPtr.GetPointer());
72  if(rval.IsNull())
73  {
74  itkExceptionMacro(<< "downcast to type " << this->GetNameOfClass()<< " failed.");
75  }
76  rval->m_Radius = m_Radius;
77  rval->m_RadiusVector = m_RadiusVector;
78  return loPtr;
79 }
80 
81 template< class TInputImage, class TOutputImage >
82 void
85 {
86  utlShowPosition(this->GetDebug());
87  if (this->m_Orientations->Rows()==0)
88  {
89  // output SH coefficients
90  utlException(this->m_Radius<0, "need to set radius");
91  utlException(this->m_SHRank<0 || this->m_RadialRank<0, "need to set this->m_SHRank and this->m_RadialRank");
92 
93  // if it is in q-space, then convert b value to q value
94  double radius = this->m_Radius;
95  if ((this->m_IsInQSpace && !this->m_IsFourier) || (!this->m_IsInQSpace && this->m_IsFourier))
96  radius = std::sqrt(this->m_Radius/(4*M_PI*M_PI*this->m_Tau));
97 
98  // std::cout << "radius = " << radius << std::endl << std::flush;
99  if (this->m_BasisType==Superclass::SPF || this->m_BasisType==Superclass::DSPF)
100  {
101  typedef SphericalPolarFourierRadialGenerator<double> SPFGenerator;
102  int n_b_sh = (this->m_SHRank+1)*(this->m_SHRank+2)/2;
103  int n_b_ra = this->m_RadialRank+1;
104  int n_b = n_b_sh * n_b_ra;
105  typename SPFGenerator::Pointer spf = SPFGenerator::New();
106 
107  // NOTE: if use ReSize, different threads modify the same data block, which is not thread safe.
108  // this->m_SPFToFeatureTransform->ReSize(n_b_sh,n_b);
109  this->m_SPFToFeatureTransform = MatrixPointer (new MatrixType(n_b_sh,n_b) );
110  this->m_SPFToFeatureTransform->Fill(0.0);
111 
112  if ( (this->m_BasisType==Superclass::SPF && !this->m_IsFourier) || (this->m_BasisType==Superclass::DSPF && this->m_IsFourier) )
113  {
114  spf->SetSPFType(SPFGenerator::SPF);
115  spf->SetScale(this->m_BasisScale);
116  for ( int n = 0; n <= this->m_RadialRank; n += 1 )
117  {
118  spf->SetN(n);
119  double spfValue = spf->Evaluate(radius);
120  int jj=0;
121  for ( int l = 0; l <= this->m_SHRank; l += 2 )
122  {
123  for ( int m = -l; m <= l; m += 1 )
124  {
125  (*this->m_SPFToFeatureTransform)(jj,n*n_b_sh+jj) = spfValue;
126  // utlPrintVar4(true, n,l,m, spf->Evaluate(radius));
127  jj++;
128  }
129  }
130  }
131  }
132  else
133  {
134  spf->SetSPFType(SPFGenerator::DSPF);
135  spf->SetScale(this->m_BasisScale);
136  for ( int n = 0; n <= this->m_RadialRank; n += 1 )
137  {
138  spf->SetN(n);
139  int jj = 0;
140  for ( int l = 0; l <= this->m_SHRank; l += 2 )
141  {
142  spf->SetL(l);
143  double spfValue = spf->Evaluate(radius);
144  for ( int m = -l; m <= l; m += 1 )
145  {
146  (*this->m_SPFToFeatureTransform)(jj,n*n_b_sh+jj) = spfValue;
147  jj++;
148  }
149  }
150  }
151  }
152  }
153  }
154  else
155  {
156  // output samples
157  if (m_RadiusVector->size()==0)
158  {
159  utlException(this->m_Radius<0, "need to set m_Radius or m_RadiusVector");
160  utl::MatchBVectorAndGradientMatrix(this->m_Radius, *m_RadiusVector, *this->m_Orientations);
161  }
162  utlException(m_RadiusVector->size()>0 && m_RadiusVector->size()!=this->m_Orientations->Rows(), "wrong size of m_RadiusVector or m_Orientations");
163  this->m_SPFIEstimator->SetBasisScale(this->m_BasisScale);
164  // std::cout << "this->m_BasisScale = " << this->m_BasisScale << std::endl << std::flush;
165  this->m_SPFIEstimator->ComputeBasisMatrix();
166  this->m_SPFToFeatureTransform = this->m_SPFIEstimator->GetBasisMatrix();
167  // this->m_SPFIEstimator->Print(std::cout<<"m_SPFIEstimator : ");
168  }
169 
170  if (this->GetDebug())
171  utl::PrintUtlMatrix(*this->m_SPFToFeatureTransform, "this->m_SPFToFeatureTransform");
172 }
173 
174 template <class TInputImage, class TOutputImage>
175 void
176 ProfileFromSPFImageFilter<TInputImage,TOutputImage>
177 ::BeforeThreadedGenerateData()
178 {
179  utlShowPosition(this->GetDebug());
180  // this->Print(std::cout<<"this=");
181  this->SetSPFIEstimator();
182 
183  if (this->m_Orientations->Rows()>0)
184  {
185  this->m_SPFIEstimator->GetSamplingSchemeQSpace()->SetOrientationsSpherical(this->m_Orientations);
186  if (m_RadiusVector->size()==0)
187  {
188  utlException(this->m_Radius<0, "need to set m_Radius or m_RadiusVector");
189  utl::MatchBVectorAndGradientMatrix(this->m_Radius, *m_RadiusVector, *this->m_Orientations);
190  }
191  this->m_SPFIEstimator->GetSamplingSchemeQSpace()->SetTau(this->m_Tau);
192  this->m_SPFIEstimator->SetSHRank(this->m_SHRank);
193  this->m_SPFIEstimator->SetRadialRank(this->m_RadialRank);
194  this->m_SPFIEstimator->SetBasisScale(this->m_BasisScale);
195  this->m_SPFIEstimator->SetScaleImage(this->m_ScaleImage);
196  this->m_SPFIEstimator->SetDebug(this->GetDebug());
197 
198  bool needToConvertToQValue = (this->m_IsInQSpace && !this->m_IsFourier) || (!this->m_IsInQSpace && this->m_IsFourier);
199  if (needToConvertToQValue)
200  this->m_SPFIEstimator->GetSamplingSchemeQSpace()->SetBVector(this->m_RadiusVector);
201  else
202  this->m_SPFIEstimator->GetSamplingSchemeQSpace()->SetRadiusVector(this->m_RadiusVector);
203  }
204 
205 
206  // this->m_SPFIEstimator->Print(std::cout<<"m_SPFIEstimator: ");
207 
208  Superclass::BeforeThreadedGenerateData();
209 }
210 
211 template <class TInputImage, class TOutputImage>
212 void
213 ProfileFromSPFImageFilter<TInputImage,TOutputImage>
214 ::ThreadedGenerateData( const typename TOutputImage::RegionType &outputRegionForThread, ThreadIdType threadId)
215 {
216  utlShowPosition(this->GetDebug());
217  typename TInputImage::ConstPointer inputPtr = this->GetInput();
218  typename TOutputImage::Pointer outputPtr = this->GetOutput();
219 
220  // Define the iterators
221  ImageRegionConstIteratorWithIndex<TInputImage> inputIt(inputPtr, outputRegionForThread);
222  ImageRegionIteratorWithIndex<TOutputImage> outputIt(outputPtr, outputRegionForThread);
223  ImageRegionIteratorWithIndex<MaskImageType> maskIt;
224  if (this->IsMaskUsed())
225  maskIt = ImageRegionIteratorWithIndex<MaskImageType>(this->m_MaskImage, outputRegionForThread);
226 
227  ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
228  typename TInputImage::PixelType inputPixel;
229  typename TInputImage::IndexType inputIndex;
230  typename TOutputImage::PixelType outputPixel;
231 
232  unsigned int outputDim = outputPtr->GetNumberOfComponentsPerPixel();
233  outputPixel.SetSize(outputDim);
234  unsigned int inputDim = inputPtr->GetNumberOfComponentsPerPixel();
235  inputPixel.SetSize(inputDim);
236 
237  inputIt.GoToBegin();
238  outputIt.GoToBegin();
239  VectorType spfVec, temp;
240 
241  Pointer selfClone = this->Clone();
242 
243  while( !inputIt.IsAtEnd() )
244  {
245  if (!this->IsMaskUsed() || (this->IsMaskUsed() && maskIt.Get()>0))
246  {
247  inputPixel = inputIt.Get();
248  if (inputPixel.GetSquaredNorm()>1e-8)
249  {
250  inputIndex = inputIt.GetIndex();
251  if (this->GetDebug())
252  std::cout << "inputIndex = " << inputIndex << std::endl << std::flush;
253  if (!IsImageEmpty(this->m_ScaleImage))
254  selfClone->SetBasisScale(this->m_ScaleImage->GetPixel(inputIndex));
255  if (selfClone->m_SPFToFeatureTransform->Size()==0)
256  selfClone->ComputeSPFToFeatureTransform();
257  spfVec = utl::VariableLengthVectorToUtlVector(inputPixel);
258  // utl::PrintUtlVector(spfVec, "spfVec");
259  // utl::PrintUtlMatrix(*selfClone->m_SPFToFeatureTransform, "*selfClone->m_SPFToFeatureTransform");
260 
261  utl::ProductUtlMv(*selfClone->m_SPFToFeatureTransform, spfVec, temp);
262  outputPixel = utl::UtlVectorToVariableLengthVector(temp);
263  if (this->GetDebug())
264  {
265  if (!IsImageEmpty(this->m_ScaleImage))
266  std::cout << "scale = " << this->m_ScaleImage->GetPixel(inputIndex) << std::endl << std::flush;
267  itk::PrintVariableLengthVector(inputPixel, "inputPixel");
268  itk::PrintVariableLengthVector(outputPixel, "outputPixel");
269  }
270  }
271  else
272  outputPixel.Fill(0.0);
273  }
274  else
275  outputPixel.Fill(0.0);
276 
277  outputIt.Set(outputPixel);
278  progress.CompletedPixel();
279  if (this->IsMaskUsed())
280  ++maskIt;
281  ++inputIt;
282  ++outputIt;
283  progress.CompletedPixel(); // potential exception thrown here
284  }
285 }
286 
287 template< class TInputImage, class TOutputImage >
288 void
289 ProfileFromSPFImageFilter< TInputImage, TOutputImage >
290 ::PrintSelf(std::ostream &os, Indent indent) const
291 {
292  Superclass::PrintSelf(os, indent);
293  os << indent
294  << ", m_Radius = " << this->m_Radius
295  << std::endl;
296  utl::PrintVector(*m_RadiusVector, "m_RadiusVector", " ", os<<indent);
297 }
298 
299 }
300 
301 
302 #endif
303 
304 
helper functions specifically used in dmritool
LightObject::Pointer InternalClone() const ITK_OVERRIDE
#define utlException(cond, expout)
Definition: utlCoreMacro.h:548
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
int RankToDimSH(const int shRank)
Definition: utlDMRI.h:193
void VerifyInputParameters() const ITK_OVERRIDE
#define utlShowPosition(cond)
Definition: utlCoreMacro.h:554
Compute some features (DWI/EAP profile, ODFs, scalar indices) from SPF coefficients.