Orfeo Toolbox  3.16
itkActiveShapeModelGradientSearchMethod.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkActiveShapeModelGradientSearchMethod.txx,v $
5  Language: C++
6  Date: $Date: 2009-04-05 23:09:19 $
7  Version: $Revision: 1.7 $
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 __itkActiveShapeModelGradientSearchMethod_txx
18 #define __itkActiveShapeModelGradientSearchMethod_txx
19 
20 #include "itkActiveShapeModelGradientSearchMethod.h"
21 
22 namespace itk
23 {
24 
25 
29 template<class TImage>
30 ActiveShapeModelGradientSearchMethod<TImage>
31 ::ActiveShapeModelGradientSearchMethod( ):m_LenghtOfProfile(3)
32 {
33  m_NumberOfIteration = 1; // default value
34  m_EigenVectors.set_size(0,0);
35  m_EigenValues.set_size(0);
36  m_MeanShape.set_size(0);
37  m_NewShape.set_size( 0);
38  m_DiffVector.set_size( 0 );
39  m_Db.set_size( 0 );
40  m_Blimit.set_size( 0 );
41  m_Valid = false;
42  m_Image = NULL;
43 }
44 
45 
49 template<class TImage>
50 void
51 ActiveShapeModelGradientSearchMethod<TImage>
52 ::GenerateData( )
53 {
54  if( !m_Image )
55  {
56  return;
57  }
58 
59  typename GradientFilterType::Pointer gradientFilter = GradientFilterType::New();
60 
61  typename Image2DType::RegionType region;
62  typename Image2DType::ConstPointer gradImage;
63 
64 
65  gradientFilter->SetInput( m_Image );
66  gradientFilter->Update();
67  gradImage = gradientFilter->GetOutput();
68 
69 
70  ConstIteratorType constIterator( gradImage, gradImage->GetLargestPossibleRegion());
71 
72  typename SampleType::Pointer sampleLandmarks = SampleType::New();
73  MeasurementVectorUIntType mv;
74  MeasurementVectorDoubleType centerShape;
75 
76  unsigned int numberOfLandmarks = m_MeanShape.size()/2;
77  unsigned int numberOfEigenValues = m_EigenValues.size();
78 
79  m_NewShape = m_MeanShape;
80  centerShape [ 0 ] = 0.0;
81  centerShape [ 1 ] = 0.0;
82 
83 
84  for ( unsigned int j = 0; j < numberOfLandmarks; j++ )
85  {
86  mv [ 0 ] = (unsigned int) (m_NewShape [ 2*j ] + 0.5);
87  mv [ 1 ] = (unsigned int) (m_NewShape [ 2*j+1 ] + 0.5);
88  sampleLandmarks->PushBack( mv );
89  centerShape [ 0 ] += mv [ 0 ];
90  centerShape [ 1 ] += mv [ 1 ];
91  }
92  centerShape [ 0 ] = centerShape [ 0 ] / numberOfLandmarks;
93  centerShape [ 1 ] = centerShape [ 1 ] / numberOfLandmarks;
94 
95  IndexType position2D, newPosition2D;
96  MeasurementVectorUIntType posRight, posLeft, dxyRef1, dxyRef2;
97  double maximumGrad, deltaX, deltaY;
98  int d, dx, dy, sx, sy;
99  unsigned int ax, ay;
100  VectorType v;
101  typename VectorType::iterator p6;
102 
103  dxyRef1.Fill(0);
104  dxyRef2.Fill(0);
105 
106  m_NewShape.set_size( 2*numberOfLandmarks );
107  m_NewShape.fill(0);
108 
109  m_DiffVector.set_size( 2*numberOfLandmarks );
110  m_DiffVector.fill(0);
111 
112  m_Db.set_size( numberOfEigenValues );
113  m_Db.fill(0);
114 
115  m_Blimit.set_size( numberOfEigenValues );
116  m_Blimit.fill(0);
117 
118 
119  for ( unsigned int i = 0; i < m_NumberOfIteration; i++ )
120  {
121  for ( unsigned int j = 0; j < numberOfLandmarks; j++ )
122  {
123  mv = sampleLandmarks->GetMeasurementVector(j);
124  position2D [ 0 ] = mv [ 0 ];
125  position2D [ 1 ] = mv [ 1 ];
126  constIterator.SetIndex(position2D);
127  maximumGrad = constIterator.Get();
128  newPosition2D = position2D;
129 
130  for (unsigned int identifier = 0; identifier<2; identifier++)
131  {
132  posRight[ identifier ] = position2D [ identifier ];
133  posLeft[ identifier ] = position2D [ identifier ];
134  }
135  if (j == 0)
136  {
137  dxyRef1 = sampleLandmarks->GetMeasurementVector(1) -
138  sampleLandmarks->GetMeasurementVector(0);
139  dxyRef2 = sampleLandmarks->GetMeasurementVector(numberOfLandmarks - 1) -
140  sampleLandmarks->GetMeasurementVector(0);
141  }
142  if (j == numberOfLandmarks - 1)
143  {
144  dxyRef1 = sampleLandmarks->GetMeasurementVector(1) -
145  sampleLandmarks->GetMeasurementVector(j);
146  dxyRef2 = sampleLandmarks->GetMeasurementVector(j-1) -
147  sampleLandmarks->GetMeasurementVector(j);
148  }
149  if ((j < numberOfLandmarks - 1) && (j > 1))
150  {
151  dxyRef1 = sampleLandmarks->GetMeasurementVector(j + 1) -
152  sampleLandmarks->GetMeasurementVector(j);
153  dxyRef2 = sampleLandmarks->GetMeasurementVector(j-1) -
154  sampleLandmarks->GetMeasurementVector(j);
155  }
156 
157  dx = dxyRef2[ 1 ] - dxyRef1[ 1 ];
158  dy = dxyRef1[ 0 ] - dxyRef2[ 0 ];
159 
160  ax = vcl_abs(dx) * 2;
161  ay = vcl_abs(dy) * 2;
162  if (dx < 0)
163  {
164  sx = -1;
165  }
166  else
167  {
168  sx = 1;
169  }
170  if (dy < 0)
171  {
172  sy = -1;
173  }
174  else
175  {
176  sy = 1;
177  }
178  unsigned int count = 0;
179  if (ax > ay)
180  {
181  d = ay - ax/2;
182  while (count < m_LenghtOfProfile)
183  {
184  if (d >= 0)
185  {
186  posRight[ 1 ] = posRight[ 1 ] + sy;
187  posLeft[ 1 ] = posLeft[ 1 ] - sy;
188  d = d - ax;
189  }
190  posRight[ 0 ] = posRight[ 0 ] + sx;
191  posLeft[ 0 ] = posLeft[ 0 ] - sx;
192  d = d + ay;
193  position2D[ 0 ] = posRight[ 0 ];
194  position2D[ 1 ] = posRight[ 1 ];
195  constIterator.SetIndex(position2D);
196  if (maximumGrad < constIterator.Get())
197  {
198  maximumGrad = constIterator.Get();
199  newPosition2D = position2D;
200  }
201  position2D[ 0 ] = posLeft[ 0 ];
202  position2D[ 1 ] = posLeft[ 1 ];
203  constIterator.SetIndex(position2D);
204  if (maximumGrad < constIterator.Get())
205  {
206  maximumGrad = constIterator.Get();
207  newPosition2D = position2D;
208  }
209  count++;
210  }
211  v.push_back (newPosition2D[ 0 ]);
212  v.push_back (newPosition2D[ 1 ]);
213  }
214  else
215  {
216  d = ax - ay/2;
217  while (count < m_LenghtOfProfile)
218  {
219  if (d >= 0)
220  {
221  posRight[ 0 ] = posRight[ 0 ] + sx;
222  posLeft[ 0 ] = posLeft[ 0 ] - sx;
223  d = d - ay;
224  }
225  posRight[ 1 ] = posRight[ 1 ] + sy;
226  posLeft[ 1 ] = posLeft[ 1 ] - sy;
227  d = d + ax;
228  position2D[ 0 ] = posRight[ 0 ];
229  position2D[ 1 ] = posRight[ 1 ];
230  constIterator.SetIndex(position2D);
231  if (maximumGrad < constIterator.Get())
232  {
233  maximumGrad = constIterator.Get();
234  newPosition2D = position2D;
235  }
236  position2D[ 0 ] = posLeft[ 0 ];
237  position2D[ 1 ] = posLeft[ 1 ];
238  constIterator.SetIndex(position2D);
239  if (maximumGrad < constIterator.Get())
240  {
241  maximumGrad = constIterator.Get();
242  newPosition2D = position2D;
243  }
244  count++;
245  }
246  v.push_back (newPosition2D[ 0 ]);
247  v.push_back (newPosition2D[ 1 ]);
248  }
249  }
250 
251  unsigned int row = 0;
252  for (p6 = v.begin(); p6 != v.end(); p6++)
253  {
254  m_NewShape[row] = (*p6);
255  row++;
256  }
257  v.erase(v.begin(), v.end());
258  sampleLandmarks->Clear();
259 
260  m_DiffVector = m_NewShape - m_MeanShape;
261 
262 
263  m_Db = m_EigenVectors.transpose() * m_DiffVector;
264 
265  for(unsigned int j = 0; j < numberOfEigenValues; j++)
266  {
267  m_Blimit[j] = 2 * vcl_sqrt(m_EigenValues[j]);
268  }
269 
270  for(unsigned int j = 0; j < numberOfEigenValues; j++)
271  {
272  if(vcl_fabs(m_Db[j]) > m_Blimit[j])
273  {
274  m_Db[j] = ( m_Db[j] / vcl_fabs(m_Db[j]) ) * m_Blimit[j];
275  }
276  }
277 
278  m_NewShape = m_MeanShape + m_EigenVectors * m_Db;
279 
280  MeasurementVectorDoubleType newCenterShape;
281  newCenterShape [ 0 ] = 0.0;
282  newCenterShape [ 1 ] = 0.0;
283 
284  for ( unsigned int j = 0; j < numberOfLandmarks; j++ )
285  {
286  mv [ 0 ] = (unsigned int) (m_NewShape [ 2*j ] + 0.5);
287  mv [ 1 ] = (unsigned int) (m_NewShape [ 2*j+1 ] + 0.5);
288  sampleLandmarks->PushBack( mv );
289  newCenterShape [ 0 ] += mv [ 0 ];
290  newCenterShape [ 1 ] += mv [ 1 ];
291  }
292  newCenterShape [ 0 ] = newCenterShape [ 0 ] / numberOfLandmarks;
293  newCenterShape [ 1 ] = newCenterShape [ 1 ] / numberOfLandmarks;
294 
295  deltaX = newCenterShape[ 0 ] - centerShape[ 0 ];
296  deltaY = newCenterShape[ 1 ] - centerShape[ 1 ];
297 
298  for ( unsigned int j = 0; j < numberOfLandmarks; j++ )
299  {
300  m_MeanShape [ 2*j ] = m_MeanShape [ 2*j ] + deltaX;
301  m_MeanShape [ 2*j+1 ] = m_MeanShape [ 2*j+1 ] + deltaY;
302  centerShape [ 0 ] += m_MeanShape [ 2*j ];
303  centerShape [ 1 ] += m_MeanShape [ 2*j+1 ];
304  }
305  centerShape [ 0 ] = centerShape [ 0 ] / numberOfLandmarks;
306  centerShape [ 1 ] = centerShape [ 1 ] / numberOfLandmarks;
307 
308  }
309 
310  /* *
311  * Remember that the moments are valid
312  */
313  m_Valid = 1;
314 
315 }
316 
317 /* *
318  * Get the new shape
319  */
320 template<class TImage>
321 typename ActiveShapeModelGradientSearchMethod<TImage>::VectorOfDoubleType
322 ActiveShapeModelGradientSearchMethod<TImage>::
323 GetNewShape()
324 {
325  if (!m_Valid)
326  {
327  throw InvalidActiveShapeModeError(__FILE__, __LINE__);
328  }
329  return m_NewShape;
330 }
331 
332 
336 template<class TImage>
337 void
338 ActiveShapeModelGradientSearchMethod<TImage>::
339 PrintSelf(std::ostream& os, Indent indent) const
340 {
341  this->Superclass::PrintSelf(os,indent);
342  os << indent << "Valid : " << m_Valid << std::endl;
343  os << indent << "Mean Shape : " << m_MeanShape << std::endl;
344  os << indent << "Eigen Vectors: " << m_EigenVectors << std::endl;
345  os << indent << "Eigen Values: " << m_EigenValues << std::endl;
346  os << indent << "Lenght of Profile: " << m_LenghtOfProfile << std::endl;
347  os << indent << "Number of Iteration: " << m_NumberOfIteration << std::endl;
348  os << indent << "Diff Vector: " << m_DiffVector << std::endl;
349  os << indent << "Db Vector: " << m_Db << std::endl;
350  os << indent << "Blimit Vector: " << m_Blimit << std::endl;
351  os << indent << "New Shape: " << m_NewShape << std::endl;
352  os << indent << "Image: " << m_Image.GetPointer() << std::endl;
353 
354  itkDebugMacro(<<" ");
355  itkDebugMacro(<<"Results of the shape model");
356  itkDebugMacro(<<"====================================");
357 
358  itkDebugMacro(<< " ");
359  itkDebugMacro(<< "================== ");
360 
361  itkDebugMacro(<< "The new shape: ");
362 
363  itkDebugMacro(<< m_NewShape);
364 
365  itkDebugMacro(<< " ");
366  itkDebugMacro(<< "+++++++++++++++++++++++++");
367 }// end PrintSelf
368 
369 } // end namespace itk
370 
371 #endif
Generated at Sat Feb 2 2013 23:23:10 for Orfeo Toolbox with doxygen 1.8.1.1