ANIMA  4.0
animaTLinksFilter.hxx
Go to the documentation of this file.
1 #pragma once
2 
3 #include "animaTLinksFilter.h"
4 
5 
6 namespace anima
7 {
8 
9 template <typename TInput, typename TOutput>
10 void TLinksFilter<TInput, TOutput>::SetInputImage(unsigned int i, const TInput* image)
11 {
12  if (i == 0)
13  this->SetNthInput(0, const_cast<TInput*>(image));
14  else
15  {
16  this->SetNthInput(m_NbInputs, const_cast<TInput*>(image));
17  m_NbInputs++;
18  }
19 
20  m_NbModalities++;
21  InConstIteratorType inputIterator (image, image->GetLargestPossibleRegion());
22  m_imagesVectorIt.push_back(inputIterator);
23  m_imagesVector.push_back(image);
24 }
25 
26 template <typename TInput, typename TOutput>
28 {
29  this->SetNthInput(m_NbInputs, const_cast<TSeedMask*>(mask));
30  m_IndexSourcesMask = m_NbInputs;
31  m_NbInputs++;
32 }
33 
34 template <typename TInput, typename TOutput>
36 {
37  this->SetNthInput(m_NbInputs, const_cast<TSeedMask*>(mask));
38  m_IndexSinksMask = m_NbInputs;
39  m_NbInputs++;
40 }
41 
42 template <typename TInput, typename TOutput>
44 {
45  this->SetNthInput(m_NbInputs, const_cast<TSeedProba*>(mask));
46  m_IndexSourcesProba = m_NbInputs;
47  m_NbInputs++;
48 }
49 
50 template <typename TInput, typename TOutput>
52 {
53  this->SetNthInput(m_NbInputs, const_cast<TSeedProba*>(mask));
54  m_IndexSinksProba = m_NbInputs;
55  m_NbInputs++;
56 }
57 
58 template <typename TInput, typename TOutput>
59 itk::Image <unsigned char,3>::ConstPointer TLinksFilter<TInput, TOutput>::GetInputSeedSourcesMask()
60 {
61  return static_cast< const TSeedMask * >
62  ( this->itk::ProcessObject::GetInput(m_IndexSourcesMask) );
63 }
64 
65 template <typename TInput, typename TOutput>
66 itk::Image <unsigned char,3>::ConstPointer TLinksFilter<TInput, TOutput>::GetInputSeedSinksMask()
67 {
68  return static_cast< const TSeedMask * >
69  ( this->itk::ProcessObject::GetInput(m_IndexSinksMask) );
70 }
71 
72 template <typename TInput, typename TOutput>
73 itk::Image <double,3>::ConstPointer TLinksFilter<TInput, TOutput>::GetInputSeedSourcesProba()
74 {
75  return static_cast< const TSeedProba * >
76  ( this->itk::ProcessObject::GetInput(m_IndexSourcesProba) );
77 }
78 
79 template <typename TInput, typename TOutput>
80 itk::Image <double,3>::ConstPointer TLinksFilter<TInput, TOutput>::GetInputSeedSinksProba()
81 {
82  return static_cast< const TSeedProba * >
83  ( this->itk::ProcessObject::GetInput(m_IndexSinksProba) );
84 }
85 
86 template <typename TInput, typename TOutput>
87 itk::DataObject::Pointer TLinksFilter<TInput, TOutput>::MakeOutput(unsigned int idx)
88 {
89  itk::DataObject::Pointer output;
90 
91  switch ( idx )
92  {
93  case 0:
94  output = ( TOutput::New() ).GetPointer();
95  break;
96  case 1:
97  output = ( TOutput::New() ).GetPointer();
98  break;
99  default:
100  std::cerr << "No output " << idx << std::endl;
101  output = NULL;
102  break;
103  }
104  return output.GetPointer();
105 }
106 
107 template <typename TInput, typename TOutput>
109 {
110  return dynamic_cast< TOutput * >( this->itk::ProcessObject::GetOutput(0) );
111 }
112 
113 template <typename TInput, typename TOutput>
115 {
116  return dynamic_cast< TOutput * >( this->itk::ProcessObject::GetOutput(1) );
117 }
118 
119 
120 template <typename TInput, typename TOutput>
121 void
124 {
125  std::cout << "Computing T-Links..." << std::endl;
126 
127  typename TOutput::Pointer output1 = this->GetOutputSources();
128  output1->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
129  output1->CopyInformation(this->GetInput(0));
130  output1->Allocate();
131  output1->FillBuffer(0);
132 
133  typename TOutput::Pointer output2 = this->GetOutputSinks();
134  output2->SetRegions(this->GetInput(0)->GetLargestPossibleRegion());
135  output2->CopyInformation(this->GetInput(0));
136  output2->Allocate();
137  output2->FillBuffer(0);
138 
139  if(this->m_TLinkMode==singleGaussianTLink && (this->GetInputSeedSourcesMask().IsNull() && this->GetInputSeedSinksMask().IsNull()))
140  {
141  std::cerr << "-- Error in T-Links computation: Single Gaussian mode requires sources or sinks mask" << std::endl;
142  exit(-1);
143  }
144  if(this->m_TLinkMode==stremTLink && (this->GetInputSeedSourcesProba().IsNull() && this->GetInputSeedSinksProba().IsNull()))
145  {
146  std::cerr << "-- Error in T-Links computation: Strem mode requires sources or sinks probabilties" << std::endl;
147  exit(-1);
148  }
149 
150  m_NbModalities = m_imagesVectorIt.size();
151 
152  switch(m_TLinkMode)
153  {
154  case singleGaussianTLink:
155  computeSingleGaussian();
156  break;
157  case stremTLink:
158  computeStrem();
159  break;
160  default:
161  return;
162  }
163 
164  OutRegionIteratorType outIterator1(output1, output1->GetLargestPossibleRegion());
165  OutRegionIteratorType outIterator2(output2, output2->GetLargestPossibleRegion());
166 
167  while (!outIterator1.IsAtEnd())
168  {
169  double val1 = m_Alpha * outIterator1.Get();
170  outIterator1.Set(val1);
171 
172  double val2 = m_Alpha * outIterator2.Get();
173  outIterator2.Set(val2);
174 
175  ++outIterator1;
176  ++outIterator2;
177  }
178 }
179 
180 template <typename TInput, typename TOutput>
181 void
184 {
185  // Just copy proba
186  std::cout << "Use of strem method..." << std::endl;
187  OutRegionIteratorType outIterator1(this->GetOutputSources(), this->GetOutputSources()->GetLargestPossibleRegion());
188  OutRegionIteratorType outIterator2(this->GetOutputSinks(), this->GetOutputSinks()->GetLargestPossibleRegion());
189 
190  SeedProbaRegionConstIteratorType sourcesIterator(this->GetInputSeedSourcesProba(), this->GetInputSeedSourcesProba()->GetLargestPossibleRegion());
191  SeedProbaRegionConstIteratorType sinksIterator(this->GetInputSeedSinksProba(), this->GetInputSeedSinksProba()->GetLargestPossibleRegion());
192 
193  while(!sourcesIterator.IsAtEnd())
194  {
195  outIterator1.Set(static_cast<OutputPixelType>(sourcesIterator.Get()));
196  outIterator2.Set(static_cast<OutputPixelType>(sinksIterator.Get()));
197 
198  ++outIterator1;
199  ++outIterator2;
200  ++sourcesIterator;
201  ++sinksIterator;
202  }
203 }
204 
205 template <typename TInput, typename TOutput>
206 void
209 {
210  if(this->GetInputSeedSourcesMask().IsNotNull())
211  {
212  std::cout << "Computing single gaussian for sources..." << std::endl;
213  computeSingleGaussianSeeds(this->GetInputSeedSourcesMask(), this->GetOutputSources(), m_MultiVarSources, this->GetInputSeedSinksMask());
214  }
215  if(this->GetInputSeedSinksMask().IsNotNull())
216  {
217  std::cout << "Computing single gaussian for sinks..." << std::endl;
218  computeSingleGaussianSeeds(this->GetInputSeedSinksMask(), this->GetOutputSinks(), m_MultiVarSinks, this->GetInputSeedSourcesMask());
219  }
220 }
221 
222 template <typename TInput, typename TOutput>
223 void
225 ::computeSingleGaussianSeeds(TSeedMask::ConstPointer seedMask, OutputImagePointer output, double multiVar, TSeedMask::ConstPointer seedMaskOpp)
226 {
227  OutRegionIteratorType outIterator(output, output->GetLargestPossibleRegion());
228  SeedMaskRegionConstIteratorType seedIterator(seedMask, seedMask->GetLargestPossibleRegion());
229 
230  double epsilon=1e-37;
231 
232  // Compute the mean and variance of each seeds for each modality
233  std::vector<double> moy, var;
234  moy.resize(m_NbModalities,0);
235  var.resize(m_NbModalities,0);
236 
237  double nbSeeds = 0.0;
238  std::vector<double> som;
239  som.resize(m_NbModalities,0);
240 
241  std::vector<double> som2;
242  som2.resize(m_NbModalities,0);
243 
244  for ( unsigned int k = 0; k < m_NbModalities; k++ )
245  {
246  m_imagesVectorIt[k].GoToBegin();
247  }
248  while(!seedIterator.IsAtEnd())
249  {
250  if(seedIterator.Get()!=0)
251  {
252  nbSeeds++;
253  for (unsigned int m = 0; m < m_NbModalities; m++)
254  {
255  som [m] += m_imagesVectorIt[m].Get();
256  som2[m] += std::pow(static_cast<double>(m_imagesVectorIt[m].Get()),2.0);
257  }
258  }
259  ++seedIterator;
260  for ( unsigned int k = 0; k < m_NbModalities; k++ )
261  {
262  ++m_imagesVectorIt[k];
263  }
264  }
265 
266  for (unsigned int m = 0; m < m_NbModalities; m++)
267  {
268  moy[m] = som[m]/nbSeeds;
269  var[m] = ((som2[m]/nbSeeds) - (moy[m]*moy[m]));
270  var[m]*=multiVar;
271  if (var[m]==0)
272  {
273  var[m]=epsilon;
274  }
275 
276  if(m_Verbose)
277  {
278  std::cout << "mean of modality " << m << ": " << moy[m] << std::endl;
279  std::cout << "variance of modality " << m << ": " << var[m] << std::endl;
280  std::cout << std::endl;
281  }
282  }
283 
284  som.clear();
285 
286  // Compute the covariances between modalities
287  std::vector<double> covar;
288  if (m_NbModalities > 1)
289  {
290  int nb_paires = 0;
291  for (unsigned int m = 0; m < m_NbModalities-1; m++)
292  {
293  for (unsigned int n = m+1; n < m_NbModalities; n++)
294  {
295  nb_paires++;
296  }
297  }
298  covar.resize(nb_paires,0);
299  som.resize(nb_paires,0);
300 
301  seedIterator.GoToBegin();
302  for ( unsigned int i = 0; i < m_NbModalities; i++ )
303  {
304  m_imagesVectorIt[i].GoToBegin();
305  }
306 
307  while(!seedIterator.IsAtEnd())
308  {
309  if(seedIterator.Get()!=0)
310  {
311  int pos=0;
312  for (unsigned int m = 0; m < m_NbModalities-1; m++)
313  {
314  for (unsigned int n = m+1; n < m_NbModalities; n++)
315  {
316  som[pos] += (m_imagesVectorIt[m].Get()-moy[m]) * (m_imagesVectorIt[n].Get()-moy[n]);
317  pos++;
318  }
319  }
320  }
321  ++seedIterator;
322  for ( unsigned int k = 0; k < m_NbModalities; k++ )
323  {
324  ++m_imagesVectorIt[k];
325 
326  }
327  }
328 
329 
330  int pos = 0;
331  for (unsigned int m = 0; m < m_NbModalities-1; m++)
332  {
333  for (unsigned int n = m+1; n < m_NbModalities; n++)
334  {
335  covar[pos] = som[pos]/(nbSeeds-1);
336  if (covar[pos]==0)
337  covar[pos]=epsilon;
338  pos++;
339  }
340  }
341  }
342 
343 
344  DoubleVariableSizeMatrixType covarMatrix;
345  covarMatrix.SetSize(m_NbModalities,m_NbModalities);
346  for (unsigned int m = 0; m < m_NbModalities; m++)
347  {
348  covarMatrix(m,m) = var[m];
349  }
350  int pos = 0;
351  for (unsigned int m = 0; m < m_NbModalities-1; m++)
352  {
353  for (unsigned int n = m+1; n < m_NbModalities; n++)
354  {
355  covarMatrix(m,n) = covarMatrix(n,m) = covar[pos];
356  pos++;
357  }
358  }
359 
360  //invert covariance matrix
361  covarMatrix=covarMatrix.GetInverse();
362 
363  // Compute the proba maps
364  MatrixTypeRes prob;
366  d.SetSize(1,m_NbModalities);
367 
368  seedIterator.GoToBegin();
369  for ( unsigned int i = 0; i < m_NbModalities; i++ )
370  {
371  m_imagesVectorIt[i].GoToBegin();
372  }
373 
374  TSeedMask::Pointer seedMaskOpp2 = TSeedMask::New();
375  seedMaskOpp2->SetRegions(seedMask->GetLargestPossibleRegion());
376  seedMaskOpp2->CopyInformation(seedMask);
377  seedMaskOpp2->Allocate();
378 
379  SeedMaskRegionIteratorType seedOppIterator2(seedMaskOpp2, seedMaskOpp2->GetLargestPossibleRegion());
380  if(seedMaskOpp.IsNull())
381  {
382  seedMaskOpp2->FillBuffer(0);
383  }
384  else
385  {
386  SeedMaskRegionConstIteratorType seedOppIterator(seedMaskOpp, seedMaskOpp->GetLargestPossibleRegion());
387  while(!seedOppIterator.IsAtEnd())
388  {
389  seedOppIterator2.Set(seedOppIterator.Get());
390  ++seedOppIterator2;
391  ++seedOppIterator;
392  }
393  }
394 
395  seedIterator.GoToBegin();
396  seedOppIterator2.GoToBegin();
397  while(!seedIterator.IsAtEnd())
398  {
399  for (unsigned int m = 0; m < m_NbModalities; m++)
400  {
401  d(0,m) = m_imagesVectorIt[m].Get() - moy[m];
402  }
403 
404  prob = d.GetVnlMatrix()*covarMatrix.GetVnlMatrix()*d.GetTranspose();
405 
406  outIterator.Set( std::exp(-0.5*prob(0,0)) );
407 
408  if(seedIterator.Get()!=0)
409  {
410  outIterator.Set(1.0);
411  }
412  if(seedOppIterator2.Get()!=0)
413  {
414  outIterator.Set(0.0);
415  }
416  ++seedOppIterator2;
417  ++outIterator;
418  ++seedIterator;
419  for ( unsigned int k = 0; k < m_NbModalities; k++ )
420  {
421  ++m_imagesVectorIt[k];
422  }
423  }
424 
425 }
426 
427 
428 } //end of namespace anima
void SetInputSeedSourcesMask(const TSeedMask *mask)
TSeedMask::ConstPointer GetInputSeedSourcesMask()
void SetInputSeedSinksMask(const TSeedMask *mask)
TSeedProba::ConstPointer GetInputSeedSinksProba()
itk::Image< PixelTypeD, 3 > TSeedProba
itk::ImageRegionConstIterator< TSeedProba > SeedProbaRegionConstIteratorType
itk::DataObject::Pointer MakeOutput(unsigned int idx)
void SetInputSeedSinksProba(const TSeedProba *mask)
TSeedMask::ConstPointer GetInputSeedSinksMask()
itk::ImageRegionConstIterator< TSeedMask > SeedMaskRegionConstIteratorType
void SetInputImage(unsigned int i, const TInput *image)
itk::ImageRegionIterator< TOutput > OutRegionIteratorType
itk::ImageRegionConstIterator< TInput > InConstIteratorType
itk::Matrix< double, 1, 1 > MatrixTypeRes
void GenerateData() ITK_OVERRIDE
itk::Image< PixelTypeUC, 3 > TSeedMask
itk::VariableSizeMatrix< NumericType > DoubleVariableSizeMatrixType
itk::ImageRegionIterator< TSeedMask > SeedMaskRegionIteratorType
TOutput::Pointer OutputImagePointer
void SetInputSeedSourcesProba(const TSeedProba *mask)
TSeedProba::ConstPointer GetInputSeedSourcesProba()
void computeSingleGaussianSeeds(TSeedMask::ConstPointer seedMask, OutputImagePointer output, double multiVar, TSeedMask::ConstPointer seedMaskOpp=ITK_NULLPTR)