Orfeo Toolbox  3.16
itkFEMImageMetricLoad.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkFEMImageMetricLoad.txx,v $
5  Language: C++
6  Date: $Date: 2009-01-29 21:28:16 $
7  Version: $Revision: 1.33 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkFEMImageMetricLoad_txx
18 #define __itkFEMImageMetricLoad_txx
19 
20 #include "itkFEMImageMetricLoad.h"
21 
22 
23 namespace itk {
24 namespace fem {
25 
26 
27 template<class TMoving,class TFixed>
28 void
29 ImageMetricLoad<TMoving , TFixed>
30 ::InitializeMetric(void)
31 {
32  if (!m_Transform) m_Transform = DefaultTransformType::New();
33  if (!m_Metric) m_Metric = DefaultMetricType::New();
34 
35  m_Temp=0.0;
36  m_Gamma=1.0;
37  m_Energy=0.0;
38 
39 //------------------------------------------------------------
40 // Set up the metric -- see MetricTest in Testing
41 //------------------------------------------------------------
42 
43  m_Metric->SetMovingImage( m_RefImage );
44  m_Metric->SetFixedImage( m_TarImage );
45 
46  typename FixedType::RegionType requestedRegion;
47  typename FixedType::SizeType size;
48  typename FixedType::IndexType tindex;
49 // typename MovingType::IndexType rindex;
50  // initialize the offset/vector part
51  for( unsigned int k = 0; k < ImageDimension; k++ )
52  {
53  //Set the size of the image region
54  size[k] = 1;
55  tindex[k]=0;
56  }
57 
58  // Set the number of integration points to zero (default for an element)
59 
60  m_NumberOfIntegrationPoints=0;
61 
62  // Set the associated region
63  requestedRegion.SetSize(size);
64  requestedRegion.SetIndex(tindex);
65  m_TarImage->SetRequestedRegion(requestedRegion);
66  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
67 
68  m_Metric->SetTransform( m_Transform.GetPointer() );
69  m_Interpolator = InterpolatorType::New();
70  m_Interpolator->SetInputImage( m_RefImage );
71  m_Metric->SetInterpolator( m_Interpolator.GetPointer() );
72 
73 
74  //------------------------------------------------------------
75  // This call is mandatory before start querying the Metric
76  // This method do all the necesary connections between the
77  // internal components: Interpolator, Transform and Images
78  //------------------------------------------------------------
79  try
80  {
81  m_Metric->Initialize();
82  }
83  catch( ExceptionObject & e )
84  {
85  std::cout << "Metric initialization failed" << std::endl;
86  std::cout << "Reason " << e.GetDescription() << std::endl;
87  }
88 
89 }
90 
91 template<class TMoving,class TFixed>
92 ImageMetricLoad<TMoving , TFixed>::ImageMetricLoad()
93 {
94  m_Metric = NULL;
95  m_Transform = NULL;
96  m_SolutionIndex = 1;
97  m_SolutionIndex2 = 0;
98  m_Sign = 1.0;
99 
100  for (unsigned int i=0; i<ImageDimension; i++)
101  {
102  m_MetricRadius[i] = 1;
103  }
104  m_MetricGradientImage = NULL;
105 
106 }
107 
108 
109 template<class TMoving,class TFixed>
110 typename ImageMetricLoad<TMoving , TFixed>::Float
111 ImageMetricLoad<TMoving , TFixed>::EvaluateMetricGivenSolution( Element::ArrayType* element,Float step)
112 {
113  Float energy=0.0,defe=0.0;
114 
115  vnl_vector_fixed<Float,2*ImageDimension> InVec(0.0);
116 
117  Element::VectorType ip,shapef;
118  Element::MatrixType solmat;
119  Element::Float w;
120 
121  Element::ArrayType::iterator elt=element->begin();
122  const unsigned int Nnodes=(*elt)->GetNumberOfNodes();
123 
124  solmat.set_size(Nnodes*ImageDimension,1);
125 
126  for(; elt!=element->end(); elt++)
127  {
128  for(unsigned int i=0; i<m_NumberOfIntegrationPoints; i++)
129  {
130  dynamic_cast<Element*>(&*(*elt))->GetIntegrationPointAndWeight(i,ip,w,m_NumberOfIntegrationPoints); // FIXME REMOVE WHEN ELEMENT NEW IS BASE CLASS
131  shapef = (*elt)->ShapeFunctions(ip);
132 
133  float solval,posval;
134  Float detJ=(*elt)->JacobianDeterminant(ip);
135 
136  for(unsigned int f=0; f<ImageDimension; f++)
137  {
138  solval=0.0;
139  posval=0.0;
140  for(unsigned int n=0; n<Nnodes; n++)
141  {
142  posval += shapef[n]*(((*elt)->GetNodeCoordinates(n))[f]);
143  float nodeval=( (m_Solution)->GetSolutionValue( (*elt)->GetNode(n)->GetDegreeOfFreedom(f) , m_SolutionIndex)
144  +(m_Solution)->GetSolutionValue( (*elt)->GetNode(n)->GetDegreeOfFreedom(f) , m_SolutionIndex2)*step);
145 
146  solval += shapef[n] * nodeval;
147  solmat[(n*ImageDimension)+f][0]=nodeval;
148  }
149  InVec[f]=posval;
150  InVec[f+ImageDimension]=solval;
151  }
152 
153  float tempe=0.0;
154  try
155  {
156  tempe=vcl_fabs(GetMetric(InVec));
157  }
158  catch( itk::ExceptionObject & )
159  {
160  // do nothing we dont care if the metric region is outside the image
161  //std::cerr << e << std::endl;
162  }
163  for(unsigned int n=0; n<Nnodes; n++)
164  {
165  itk::fem::Element::Float temp=shapef[n]*tempe*w*detJ;
166  energy += temp;
167  }
168  }
169 
170  defe += 0.0;//(double)(*elt)->GetElementDeformationEnergy( solmat );
171  }
172 
173  //std::cout << " def e " << defe << " sim e " << energy*m_Gamma << std::endl;
174  return vcl_fabs((double)energy*(double)m_Gamma-(double)defe);
175 
176 }
177 
178 template<class TMoving,class TFixed>
179 typename ImageMetricLoad<TMoving , TFixed>::VectorType
180 ImageMetricLoad<TMoving , TFixed>::Fe
181 ( VectorType Gpos, VectorType Gsol)
182 {
183  // We assume the vector input is of size 2*ImageDimension.
184  // The 0 to ImageDimension-1 elements contain the position, p,
185  // in the reference image. The next ImageDimension to 2*ImageDimension-1
186  // elements contain the value of the vector field at that point, v(p).
187  //
188  // Thus, we evaluate the derivative at the point p+v(p) with respect to
189  // some region of the target (fixed) image by calling the metric with
190  // the translation parameters as provided by the vector field at p.
191  //------------------------------------------------------------
192  // Set up transform parameters
193  //------------------------------------------------------------
194 
195  VectorType OutVec;
196 
197  for( unsigned int k = 0; k < ImageDimension; k++ )
198  {
199  if ( vnl_math_isnan(Gpos[k]) || vnl_math_isinf(Gpos[k]) ||
200  vnl_math_isnan(Gsol[k]) || vnl_math_isinf(Gsol[k]) ||
201  vcl_fabs(Gpos[k]) > 1.e33 || vcl_fabs(Gsol[k]) > 1.e33 )
202  {
203  OutVec.set_size(ImageDimension); OutVec.fill(0.0); return OutVec;
204  }
205  }
206 
207  ParametersType parameters( m_Transform->GetNumberOfParameters() );
208  typename FixedType::RegionType requestedRegion;
209  FixedRadiusType regionRadius;
210  typename FixedType::IndexType tindex;
211  typename MovingType::IndexType rindex;
212  OutVec.set_size(ImageDimension);
213 
214  int lobordercheck=0,hibordercheck=0;
215  for( unsigned int k = 0; k < ImageDimension; k++ )
216  {
217  //Set the size of the image region
218  parameters[k]= Gsol[k]; // this gives the translation by the vector field
219  rindex[k] =(long)(Gpos[k]+Gsol[k]+0.5); // where the piece of reference image currently lines up under the above translation
220  tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2; // position in reference image
221  hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
222  lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
223  if (hibordercheck >= 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
224  else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
225  else regionRadius[k]=m_MetricRadius[k];
226  tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2; // position in reference image
227  }
228 
229  // Set the associated region
230 
231  requestedRegion.SetSize(regionRadius);
232  requestedRegion.SetIndex(tindex);
233 
234  m_TarImage->SetRequestedRegion(requestedRegion);
235  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
236 
237  //--------------------------------------------------------
238  // Get metric values
239 
240  typename MetricBaseType::MeasureType measure;
241  typename MetricBaseType::DerivativeType derivative;
242 
243  try
244  {
245  m_Metric->GetValueAndDerivative( parameters, measure, derivative );
246  // m_Metric->GetDerivative( parameters, derivative );
247  }
248  catch( ... )
249  {
250  // do nothing we don't care if the metric lies outside the image sometimes
251  //std::cerr << e << std::endl;
252  }
253 
254  m_Energy += (double)measure;
255  float gmag=0.0;
256  for( unsigned int k = 0; k < ImageDimension; k++ )
257  {
258  if (lobordercheck < 0 || hibordercheck >=0 ||
259  vnl_math_isnan(derivative[k]) || vnl_math_isinf(derivative[k]) )
260  {
261  OutVec[k]=0.0;
262  }
263  else OutVec[k]= m_Sign*m_Gamma*derivative[k];
264  gmag += OutVec[k]*OutVec[k];
265  }
266  if (gmag==0.0) gmag=1.0;
267  // NOTE : POSSIBLE THAT DERIVATIVE DIRECTION POINTS UP OR DOWN HILL!
268  // IN FACT, IT SEEMS MEANSQRS AND NCC POINT IN DIFFT DIRS
269  //std::cout << " deriv " << derivative << " val " << measure << endl;
270  //if (m_Temp !=0.0)
271  //return OutVec * vcl_exp(-1.*OutVec.magnitude()/m_Temp);
272  //else
273  return OutVec/vcl_sqrt(gmag);
274 }
275 
276 
277 template<class TMoving,class TFixed>
278 typename ImageMetricLoad<TMoving , TFixed>::Float
279 ImageMetricLoad<TMoving , TFixed>::GetMetric
280 ( VectorType InVec)
281 {
282  // We assume the vector input is of size 2*ImageDimension.
283  // The 0 to ImageDimension-1 elements contain the position, p,
284  // in the reference image. The next ImageDimension to 2*ImageDimension-1
285  // elements contain the value of the vector field at that point, v(p).
286  //
287  // Thus, we evaluate the derivative at the point p+v(p) with respect to
288  // some region of the target (fixed) image by calling the metric with
289  // the translation parameters as provided by the vector field at p.
290  //------------------------------------------------------------
291  // Set up transform parameters
292  //------------------------------------------------------------
293  ParametersType parameters( m_Transform->GetNumberOfParameters());
294  typename FixedType::RegionType requestedRegion;
295  typename FixedType::IndexType tindex;
296  typename MovingType::IndexType rindex;
297  FixedRadiusType regionRadius;
298  VectorType OutVec(ImageDimension,0.0); // gradient direction
299  //std::cout << " pos translation " << InVec << endl;
300  // initialize the offset/vector part
301  for( unsigned int k = 0; k < ImageDimension; k++ )
302  {
303  //Set the size of the image region
304  parameters[k]= InVec[k+ImageDimension]; // this gives the translation by the vector field
305  rindex[k] =(long)(InVec[k]+InVec[k+ImageDimension]+0.5); // where the piece of reference image currently lines up under the above translation
306  tindex[k]= (long)(InVec[k]+0.5)-(long)m_MetricRadius[k]/2; // position in reference image
307  int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
308  int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
309  if (hibordercheck > 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
310  else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
311  else regionRadius[k]=m_MetricRadius[k];
312  tindex[k]= (long)(InVec[k]+0.5)-(long)regionRadius[k]/2; // position in reference image
313  }
314 
315  // Set the associated region
316 
317  requestedRegion.SetSize(regionRadius);
318  requestedRegion.SetIndex(tindex);
319 
320  m_TarImage->SetRequestedRegion(requestedRegion);
321  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
322 
323  //--------------------------------------------------------
324  // Get metric values
325 
326  typename MetricBaseType::MeasureType measure=0.0;
327  try
328  {
329  measure=m_Metric->GetValue( parameters);
330  }
331  catch( ... )
332  {
333  // do nothing we dont care if the metric lies outside the image sometimes
334  //std::cerr << e << std::endl;
335  }
336 
337 
338  return (Float) measure;
339 }
340 
341 
342 template<class TMoving,class TFixed>
343 typename ImageMetricLoad<TMoving , TFixed>::VectorType
344 ImageMetricLoad<TMoving , TFixed>::MetricFiniteDiff
345 ( VectorType Gpos,
346  VectorType Gsol )
347 {
348 
349  typename MetricBaseType::MeasureType measure;
350  ParametersType parameters( ImageDimension );
351  typename FixedType::RegionType requestedRegion;
352  typename FixedType::IndexType tindex;
353  FixedRadiusType regionRadius;
354 
355  VectorType OutVec;
356  OutVec.set_size(ImageDimension);
357 
358  for( unsigned int k = 0; k < ImageDimension; k++ )
359  {
360  parameters[k]= Gsol[k]; // this gives the translation by the vector field
361  tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2; // position in reference image
362  if (tindex[k] > m_TarSize[k]-1 || tindex[k] < 0) tindex[k]=(long)(Gpos[k]+0.5);
363  int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
364  int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
365  if (hibordercheck >= 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
366  else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
367  else regionRadius[k]=m_MetricRadius[k];
368  tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2; // position in reference image
369  }
370 
371  unsigned int row;
372  typename ImageType::IndexType difIndex[ImageDimension][2];
373 
374  typename MetricBaseType::MeasureType dPixL,dPixR;
375  for(row=0; row< ImageDimension;row++)
376  {
377  difIndex[row][0]=tindex;
378  difIndex[row][1]=tindex;
379  if (tindex[row] < m_TarSize[row]-1)
380  {
381  difIndex[row][0][row]=tindex[row]+1;
382  }
383  if (tindex[row] > 0 )
384  {
385  difIndex[row][1][row]=tindex[row]-1;
386  }
387  try
388  {
389  requestedRegion.SetIndex(difIndex[row][1]);
390  requestedRegion.SetSize(regionRadius);
391  m_TarImage->SetRequestedRegion(requestedRegion);
392  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
393  dPixL=m_Metric->GetValue( parameters);
394  }
395  catch( ... )
396  {
397  dPixL=0.0;
398  }
399  try
400  {
401  requestedRegion.SetIndex(difIndex[row][0]);
402  requestedRegion.SetSize(regionRadius);
403  m_TarImage->SetRequestedRegion(requestedRegion);
404  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
405  dPixR=m_Metric->GetValue( parameters);
406  }
407  catch( ... )
408  {
409  dPixR=0.0;
410  }
411 
412  OutVec[row]=dPixL-dPixR;
413  }
414  return OutVec;
415 }
416 
417 
418 template<class TMoving,class TFixed>
419 typename ImageMetricLoad<TMoving , TFixed>::VectorType
420 ImageMetricLoad<TMoving , TFixed>::GetPolynomialFitToMetric
421 ( VectorType Gpos,
422  VectorType Gsol )
423 {
424 
425  //discrete orthogonal polynomial fitting
426  //see p.394-403 haralick computer and robot vision
427  //
428  //here, use chebyshev polynomials for fitting a plane to the data
429  //
430  //f(x,y,z) = a0 + a1*x + a2*y + a3*z
431  //
432  ParametersType parameters( ImageDimension );
433  typename FixedType::RegionType requestedRegion;
434  typename FixedType::IndexType tindex;
435  FixedRadiusType regionRadius;
436 
437  typename ImageType::IndexType temp;
438 
439  VectorType chebycoefs; // gradient direction
440  chebycoefs.set_size(ImageDimension);
441  double chebycoefs0=0.0; // the constant term
442  double datatotal=0.0;
443  double a0norm=1.0;
444  double a1norm=1.0/2.0;
445 
446  double met, ind1,ind2;
447  double inds[3]; inds[0]=-1.0; inds[1]=0.0; inds[2]=1.0;
448 
449  for( unsigned int k = 0; k < ImageDimension; k++ )
450  {
451  a0norm /= 3.0;
452  if (k < ImageDimension-1)
453  {
454  a1norm /= 3.0;
455  }
456  chebycoefs[k]=0.0;
457  parameters[k]= Gsol[k]; // this gives the translation by the vector field
458  tindex[k]= (long)(Gpos[k]+0.5)-(long)m_MetricRadius[k]/2; // position in reference image
459  if (tindex[k] > m_TarSize[k]-1 || tindex[k] < 0) tindex[k]=(long)(Gpos[k]+0.5);
460  int hibordercheck=(int)tindex[k]+(int)m_MetricRadius[k]-(int)m_TarSize[k];
461  int lobordercheck=(int)tindex[k]-(int)m_MetricRadius[k];
462  if (hibordercheck >= 0) regionRadius[k]=m_MetricRadius[k]-(long)hibordercheck-1;
463  else if (lobordercheck < 0) regionRadius[k]=m_MetricRadius[k]+(long)lobordercheck;
464  else regionRadius[k]=m_MetricRadius[k];
465  tindex[k]= (long)(Gpos[k]+0.5)-(long)regionRadius[k]/2; // position in reference image
466  }
467 
468 
469  if (ImageDimension==2)
470  {
471 
472  double measure[3][3];
473  for(int row=-1; row< 2; row++)
474  {
475  for(int col=-1; col< 2; col++)
476  {
477 
478  temp[0]=tindex[0]+(long)row;
479  temp[1]=tindex[1]+(long)col;
480 
481  for (unsigned int i=0; i<ImageDimension; i++)
482  {
483  if (temp[i] > m_TarSize[i]-1)
484  {
485  temp[i]=m_TarSize[i]-1;
486  }
487  else if (temp[i] < 0 )
488  {
489  temp[i]=0;
490  }
491  }
492 
493  requestedRegion.SetIndex(temp);
494  requestedRegion.SetSize(regionRadius);
495  m_TarImage->SetRequestedRegion(requestedRegion);
496  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
497  measure[row+1][col+1]=0.0;
498 
499  try
500  {
501  measure[row+1][col+1]=m_Metric->GetValue( parameters);
502  }
503  catch( ... )
504  {
505  }
506 
507 
508  datatotal += measure[row+1][col+1];
509  }
510  }
511  for( unsigned int cb1 = 0; cb1 < 3; cb1++ )
512  {
513  for( unsigned int cb2 = 0; cb2 < 3; cb2++ )
514  {
515  met=measure[cb1][cb2];
516  ind1=inds[cb1]*a1norm;
517  ind2=inds[cb2]*a1norm;
518  chebycoefs[0] += met*ind1;
519  chebycoefs[1] += met*ind2;
520  }
521  }
522  }
523  else if (ImageDimension == 3)
524  {
525 
526  double measure3D[3][3][3];
527  for(int row=-1; row< 2; row++)
528  {
529  for(int col=-1; col< 2; col++)
530  {
531  for(int z=-1; z< 2; z++)
532  {
533 
534  temp[0]=tindex[0]+(long)row;
535  temp[1]=tindex[1]+(long)col;
536  temp[2]=tindex[2]+(long)z;
537 
538  for (unsigned int i=0; i<ImageDimension; i++)
539  {
540  if (temp[i] > m_TarSize[i]-1)
541  {
542  temp[i]=m_TarSize[i]-1;
543  }
544  else if (temp[i] < 0 )
545  {
546  temp[i]=0;
547  }
548  }
549 
550  requestedRegion.SetIndex(temp);
551  requestedRegion.SetSize(regionRadius);
552  m_TarImage->SetRequestedRegion(requestedRegion);
553  m_Metric->SetFixedImageRegion( m_TarImage->GetRequestedRegion() );
554  measure3D[row+1][col+1][z+1]=0.0;
555 
556  try
557  {
558  measure3D[row+1][col+1][z+1]=m_Metric->GetValue( parameters);
559  }
560  catch( ... )
561  {
562  }
563 
564 
565  datatotal += measure3D[row+1][col+1][z+1];
566  }
567  }
568  }
569  for( unsigned int cb1 = 0; cb1 < 2; cb1++ )
570  {
571  for( unsigned int cb2 = 0; cb2 < 2; cb2++ )
572  {
573  for( unsigned int cb3 = 0; cb3 < 2; cb3++ )
574  {
575  chebycoefs[0] += measure3D[cb1][cb2][cb3]*inds[cb1]*a1norm;
576  chebycoefs[1] += measure3D[cb1][cb2][cb3]*inds[cb2]*a1norm;
577  chebycoefs[2] += measure3D[cb1][cb2][cb3]*inds[cb3]*a1norm;
578  }
579  }
580  }
581  }
582 
583  chebycoefs0=a0norm*datatotal;
584 // std::cout << " cb " << chebycoefs << std::endl;
585  return chebycoefs;
586 
587 }
588 
589 
590 template<class TMoving,class TFixed>
591 int ImageMetricLoad<TMoving,TFixed>::CLID()
592 {
593  static const int CLID_ = FEMOF::Register( ImageMetricLoad::NewB,(std::string("ImageMetricLoad(")
594  +typeid(TMoving).name()+","+typeid(TFixed).name()+")").c_str());
595  return CLID_;
596 }
597 
598 
599 template<class TMoving,class TFixed>
600 const int ImageMetricLoad<TMoving,TFixed>::m_DummyCLID=ImageMetricLoad<TMoving,TFixed>::CLID();
601 
602 
603 } // end namespace fem
604 } // end namespace itk
605 
606 #endif
Generated at Sat Feb 2 2013 23:36:46 for Orfeo Toolbox with doxygen 1.8.1.1