Orfeo Toolbox  3.16
itkDeformableMesh3DFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkDeformableMesh3DFilter.txx,v $
5  Language: C++
6  Date: $Date: 2009-09-17 11:04:01 $
7  Version: $Revision: 1.65 $
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 __itkDeformableMesh3DFilter_txx
18 #define __itkDeformableMesh3DFilter_txx
19 
20 #ifdef _MSC_VER
21 #pragma warning ( disable : 4786 )
22 #endif
23 
24 #include "itkDeformableMesh3DFilter.h"
25 #include "itkNumericTraits.h"
26 
27 namespace itk
28 {
29 /* Constructor. */
30 template <typename TInputMesh, typename TOutputMesh>
31 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
32 ::DeformableMesh3DFilter()
33 {
34  m_Step = 0;
35  m_StepThreshold = 0;
36  m_PotentialOn = 0;
37  m_K = 0;
38  m_ObjectLabel = 0;
39  m_Scale.Fill( 1.0 );
40  m_Stiffness.Fill( 0.1 );
41  m_TimeStep = 0.01;
42  m_PotentialMagnitude = NumericTraits<PixelType>::One;
43  m_GradientMagnitude = NumericTraits<PixelType>::One;
44 
45  m_ImageDepth = 0;
46  m_ImageHeight = 0;
47  m_ImageWidth = 0;
48 
49  typename TOutputMesh::Pointer output = TOutputMesh::New();
50  this->ProcessObject::SetNumberOfRequiredOutputs(1);
51  this->ProcessObject::SetNthOutput(0, output.GetPointer());
52 }
53 
54 template <typename TInputMesh, typename TOutputMesh>
55 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
56 ::~DeformableMesh3DFilter()
57 {
58  if (m_K)
59  {
60  delete [] (m_K);
61  m_K = 0;
62  }
63 }
64 
65 
66 /* PrintSelf. */
67 template <typename TInputMesh, typename TOutputMesh>
68 void
69 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
70 ::PrintSelf(std::ostream& os, Indent indent) const
71 {
72  Superclass::PrintSelf(os,indent);
73  os << indent << "Stiffness = " << m_Stiffness;
74  os << indent << "PotentialMagnitude = "
75  << static_cast<typename NumericTraits<PixelType>::PrintType>(m_PotentialMagnitude)
76  << std::endl;
77  os << indent << "GradientMagnitude = "
78  << static_cast<typename NumericTraits<PixelType>::PrintType>(m_GradientMagnitude)
79  << std::endl;
80  os << indent << "Scale = " << m_Scale;
81  os << indent << "TimeStep = " << m_TimeStep << std::endl;
82  os << indent << "PotentialOn = " << m_PotentialOn << std::endl;
83  os << indent << "ObjectLabel = "
84  << static_cast<typename NumericTraits<unsigned char>::PrintType>(m_ObjectLabel)
85  << std::endl;
86  os << indent << "StepThreshold = " << m_StepThreshold << std::endl;
87  if (m_Normals)
88  {
89  os << indent << "Normals = " << m_Normals << std::endl;
90  }
91  else
92  {
93  os << indent << "Normals = " << "(None)" << std::endl;
94  }
95  if (m_Gradient)
96  {
97  os << indent << "Gradient = " << m_Gradient << std::endl;
98  }
99  else
100  {
101  os << indent << "Gradient = " << "(None)" << std::endl;
102  }
103  os << indent << "Step = " << m_Step << std::endl;
104  os << indent << "ImageDepth = " << m_ImageDepth << std::endl;
105  os << indent << "ImageHeight = " << m_ImageHeight << std::endl;
106  os << indent << "ImageWidth = " << m_ImageWidth << std::endl;
107 }/* End PrintSelf. */
108 
109 /* Set default value of parameters and initialize local data container
110  * such as forces, displacements and displacement derivatives. */
111 template <typename TInputMesh, typename TOutputMesh>
112 void
113 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
114 ::Initialize()
115 {
116  m_NumberOfNodes = this->GetInput(0)->GetNumberOfPoints();
117  m_NumberOfCells = this->GetInput(0)->GetNumberOfCells();
118 
119  m_Forces = InputMeshType::New();
120  m_Displacements = InputMeshType::New();
121  m_Derives = InputMeshType::New();
122  m_Normals = InputMeshType::New();
123  m_Locations = InputMeshType::New();
124 
125  InputMeshConstPointer inputMesh = this->GetInput(0);
126  InputPointsContainerConstPointer myPoints = inputMesh->GetPoints();
127  InputPointsContainerConstIterator points = myPoints->Begin();
128 
129  InputPointsContainerPointer myForces = m_Forces->GetPoints();
130  myForces->Reserve(m_NumberOfNodes);
131  InputPointsContainerIterator forces = myForces->Begin();
132 
133  InputPointsContainerPointer myDerives = m_Derives->GetPoints();
134  myDerives->Reserve(m_NumberOfNodes);
135  InputPointsContainerIterator derives = myDerives->Begin();
136 
137  InputPointsContainerPointer myDisplacements = m_Displacements->GetPoints();
138  myDisplacements->Reserve(m_NumberOfNodes);
139  InputPointsContainerIterator displacements = myDisplacements->Begin();
140 
141  InputPointsContainerPointer myNormals = m_Normals->GetPoints();
142  myNormals->Reserve(m_NumberOfNodes);
143  InputPointsContainerIterator normals = myNormals->Begin();
144 
145  InputPointsContainerPointer myLocations = m_Locations->GetPoints();
146  myLocations->Reserve(m_NumberOfNodes);
147  InputPointsContainerIterator locations = myLocations->Begin();
148 
149  InputCellsContainerConstPointer myCells = inputMesh->GetCells();
150  InputCellsContainerConstIterator cells = myCells->Begin();
151 
152  ImageSizeType ImageSize = m_Gradient->GetBufferedRegion().GetSize();
153 
154  //---------------------------------------------------------------------
155  //Get the image width/height and depth
156  //---------------------------------------------------------------------
157  m_ImageWidth = ImageSize[0];
158  m_ImageHeight = ImageSize[1];
159  m_ImageDepth = ImageSize[2];
160 
161  m_ModelXUpLimit = 0;
162  m_ModelXDownLimit = m_ImageWidth;
163  m_ModelYUpLimit = 0;
164  m_ModelYDownLimit = m_ImageHeight;
165  m_ModelZUpLimit = 0;
166  m_ModelZDownLimit = m_ImageDepth;
167 
168  InputPointType d; d.Fill(0.0);
169 
170  int j = 0;
171  while( points != myPoints->End() )
172  {
173  for (int i=0; i<3; i++ )
174  {
175  locations.Value()[i] = m_Scale[i] * points.Value()[i];
176  }
177  ++points;
178  ++locations;
179  forces.Value() = d;
180  ++forces;
181  normals.Value() = d;
182  ++normals;
183  derives.Value() = d;
184  ++derives;
185  displacements.Value() = d;
186  ++displacements;
187  m_Forces->SetPointData(j, 0.0);
188  j++;
189  }
190 
191  const unsigned long *tp;
192  PixelType x = 0.0;
193  PixelType* x_pt;
194  x_pt = &x;
195  while( cells != myCells->End() )
196  {
197  typename InputMeshType::CellType * cellPtr = cells.Value();
198  tp = cellPtr->GetPointIds();
199  for ( int i=0; i<3; i++ ) {
200  m_Forces->GetPointData((int)(tp[i]), x_pt);
201  x = x + 1.0;
202  m_Forces->SetPointData((int)(tp[i]), x);
203  }
204  ++cells;
205  }
206 
207  this->SetDefaultStiffnessMatrix();
208 
209  // This prevents unnecessary re-executions of the pipeline.
210  OutputMeshPointer outputMesh = this->GetOutput();
211  outputMesh->SetBufferedRegion( outputMesh->GetRequestedRegion() );
212 }
213 
214 /* Set the stiffness matrix. */
215 template <typename TInputMesh, typename TOutputMesh>
216 void
217 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
218 ::SetDefaultStiffnessMatrix ()
219 {
220  double us = 0.5;
221  double vs = 0.5;
222  double a = us*us, b = vs*vs;
223  double area = us*vs/2, k00, k01, k02, k11, k12, k22;
224 
225  k00 = area * (m_Stiffness[1]/a + m_Stiffness[1]/b + m_Stiffness[0]);
226  k01 = area * (-m_Stiffness[1]/a + m_Stiffness[0] * 0.5);
227  k11 = area * (m_Stiffness[1]/a + m_Stiffness[0]);
228  k02 = area * (-m_Stiffness[1]/b + m_Stiffness[0] * 0.5);
229  k12 = area * m_Stiffness[0] * 0.5;
230  k22 = area * (m_Stiffness[1]/b + m_Stiffness[0]);
231 
232  m_StiffnessMatrix[0][0][0] = k00;
233  m_StiffnessMatrix[0][0][1] = k01;
234  m_StiffnessMatrix[0][0][2] = k02;
235  m_StiffnessMatrix[0][0][3] = 0.0;
236  m_StiffnessMatrix[0][1][0] = k01;
237  m_StiffnessMatrix[0][1][1] = k11;
238  m_StiffnessMatrix[0][1][2] = k12;
239  m_StiffnessMatrix[0][1][3] = 0.0;
240  m_StiffnessMatrix[0][2][0] = k02;
241  m_StiffnessMatrix[0][2][1] = k12;
242  m_StiffnessMatrix[0][2][2] = k22;
243  m_StiffnessMatrix[0][2][3] = 0.0;
244  m_StiffnessMatrix[0][3][0] = 0.0;
245  m_StiffnessMatrix[0][3][1] = 0.0;
246  m_StiffnessMatrix[0][3][2] = 0.0;
247  m_StiffnessMatrix[0][3][3] = 1.0;
248 
249 }
250 
252 template <typename TInputMesh, typename TOutputMesh>
253 void
254 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
255 ::SetStiffnessMatrix ( StiffnessMatrixType *stiff, int i )
256 {
257  m_StiffnessMatrix[i][0][0] = stiff[0][0];
258  m_StiffnessMatrix[i][0][1] = stiff[0][1];
259  m_StiffnessMatrix[i][0][2] = stiff[0][2];
260  m_StiffnessMatrix[i][0][3] = stiff[0][3];
261  m_StiffnessMatrix[i][1][0] = stiff[1][0];
262  m_StiffnessMatrix[i][1][1] = stiff[1][1];
263  m_StiffnessMatrix[i][1][2] = stiff[1][2];
264  m_StiffnessMatrix[i][1][3] = stiff[1][3];
265  m_StiffnessMatrix[i][2][0] = stiff[2][0];
266  m_StiffnessMatrix[i][2][1] = stiff[2][1];
267  m_StiffnessMatrix[i][2][2] = stiff[2][2];
268  m_StiffnessMatrix[i][2][3] = stiff[2][3];
269  m_StiffnessMatrix[i][3][0] = stiff[3][0];
270  m_StiffnessMatrix[i][3][1] = stiff[3][1];
271  m_StiffnessMatrix[i][3][2] = stiff[3][2];
272  m_StiffnessMatrix[i][3][3] = stiff[3][3];
273 }
274 
276 template <typename TInputMesh, typename TOutputMesh>
277 void
278 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
279 ::SetMeshStiffness ()
280 {
281  InputMeshConstPointer inputMesh = this->GetInput(0);
282  InputCellDataContainerConstPointer myCellData = inputMesh->GetCellData();
283  InputCellDataContainerConstIterator celldata = myCellData->Begin();
284 
285  m_K = new StiffnessMatrixRawPointer[m_NumberOfCells];
286 
287  unsigned int j = 0;
288  while (celldata != myCellData->End())
289  {
290  const double x = celldata.Value();
291  m_K[j] = m_StiffnessMatrix+((int) x);
292  ++celldata;
293  j++;
294  }
295 }
296 
297 
298 /* Compute the derivatives using d' + Kd = f. */
299 template <typename TInputMesh, typename TOutputMesh>
300 void
301 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
302 ::ComputeDt()
303 {
304  int i;
305  const unsigned long *tp;
306 
307  InputMeshConstPointer inputMesh = this->GetInput(0);
308  InputCellsContainerConstPointer myCells = inputMesh->GetCells();
309  InputCellsContainerConstIterator cells = myCells->Begin();
310 
311  InputPointsContainerPointer myForces = m_Forces->GetPoints();
312  InputPointsContainerIterator forces = myForces->Begin();
313 
314  InputPointsContainerPointer myDerives = m_Derives->GetPoints();
315  InputPointsContainerIterator derives = myDerives->Begin();
316 
317  double p = 1.0;
318  i = 0;
319  InputPointType v1, v2, v3;
320  v1.Fill(0.);
321  v2.Fill(0.);
322  v3.Fill(0.);
323 
324  while( cells != myCells->End() )
325  {
326  tp = cells.Value()->GetPointIds();
327  ++cells;
328  m_Displacements->GetPoint (tp[0], &v1);
329  m_Displacements->GetPoint (tp[1], &v2);
330  m_Displacements->GetPoint (tp[2], &v3);
331  v1[0] *= m_K[i]->get(0, 0)*p;
332  v1[1] *= m_K[i]->get(0, 0)*p;
333  v1[2] *= m_K[i]->get(0, 0)*p;
334  v2[0] *= m_K[i]->get(0, 1)*p;
335  v2[1] *= m_K[i]->get(0, 1)*p;
336  v2[2] *= m_K[i]->get(0, 1)*p;
337  v3[0] *= m_K[i]->get(0, 2)*p;
338  v3[1] *= m_K[i]->get(0, 2)*p;
339  v3[2] *= m_K[i]->get(0, 2)*p;
340  v1[0] += v2[0]+v3[0];
341  v1[1] += v2[1]+v3[1];
342  v1[2] += v2[2]+v3[2];
343  m_Forces->GetPoint (tp[0], &v2);
344 
345  v2[0] -= v1[0];
346  v2[1] -= v1[1];
347  v2[2] -= v1[2];
348 
349  m_Forces->SetPoint (tp[0], v2);
350 
351  m_Displacements->GetPoint (tp[0], &v1);
352  m_Displacements->GetPoint (tp[1], &v2);
353  m_Displacements->GetPoint (tp[2], &v3);
354  v1[0] *= m_K[i]->get(1, 0)*p;
355  v1[1] *= m_K[i]->get(1, 0)*p;
356  v1[2] *= m_K[i]->get(1, 0)*p;
357  v2[0] *= m_K[i]->get(1, 1)*p;
358  v2[1] *= m_K[i]->get(1, 1)*p;
359  v2[2] *= m_K[i]->get(1, 1)*p;
360  v3[0] *= m_K[i]->get(1, 2)*p;
361  v3[1] *= m_K[i]->get(1, 2)*p;
362  v3[2] *= m_K[i]->get(1, 2)*p;
363  v1[0] += v2[0]+v3[0];
364  v1[1] += v2[1]+v3[1];
365  v1[2] += v2[2]+v3[2];
366  m_Forces->GetPoint (tp[1], &v2);
367 
368  v2[0] -= v1[0];
369  v2[1] -= v1[1];
370  v2[2] -= v1[2];
371 
372  m_Forces->SetPoint (tp[1], v2);
373 
374  m_Displacements->GetPoint (tp[0], &v1);
375  m_Displacements->GetPoint (tp[1], &v2);
376  m_Displacements->GetPoint (tp[2], &v3);
377  v1[0] *= m_K[i]->get(2, 0)*p;
378  v1[1] *= m_K[i]->get(2, 0)*p;
379  v1[2] *= m_K[i]->get(2, 0)*p;
380  v2[0] *= m_K[i]->get(2, 1)*p;
381  v2[1] *= m_K[i]->get(2, 1)*p;
382  v2[2] *= m_K[i]->get(2, 1)*p;
383  v3[0] *= m_K[i]->get(2, 2)*p;
384  v3[1] *= m_K[i]->get(2, 2)*p;
385  v3[2] *= m_K[i]->get(2, 2)*p;
386  v1[0] += v2[0]+v3[0];
387  v1[1] += v2[1]+v3[1];
388  v1[2] += v2[2]+v3[2];
389  m_Forces->GetPoint (tp[2], &v2);
390 
391  v2[0] -= v1[0];
392  v2[1] -= v1[1];
393  v2[2] -= v1[2];
394 
395  m_Forces->SetPoint (tp[2], v2);
396  ++i;
397  }
398 
399  while ( derives != myDerives->End() )
400  {
401  derives.Value() = forces.Value();
402  ++derives;
403  ++forces;
404  }
405 }
406 
408 template <typename TInputMesh, typename TOutputMesh>
409 void
410 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
411 ::Advance()
412 {
413  typename TInputMesh::PointType s, d, ds;
414 
415  InputPointsContainerPointer myDerives = m_Derives->GetPoints();
416  InputPointsContainerIterator derives = myDerives->Begin();
417  InputPointsContainerPointer myDisplacements = m_Displacements->GetPoints();
418  InputPointsContainerIterator displacements = myDisplacements->Begin();
419  InputPointsContainerPointer myPoints = m_Locations->GetPoints();
420  InputPointsContainerIterator points = myPoints->Begin();
421 
422  while( derives != myDerives->End() )
423  {
424  ds = derives.Value();
425  s = points.Value();
426  d = displacements.Value();
427  s[0] += m_TimeStep*ds[0];
428  s[1] += m_TimeStep*ds[1];
429  s[2] += m_TimeStep*ds[2];
430  d[0] += m_TimeStep*ds[0];
431  d[1] += m_TimeStep*ds[1];
432  d[2] += m_TimeStep*ds[2];
433 
435  if ( (s[0] > 0) && (s[1] > 0) && (s[2] > 0) &&
436  (s[2] < m_ImageDepth) && (s[0] < m_ImageWidth) && (s[1] < m_ImageHeight) )
437  {
438  points.Value() = s;
439  displacements.Value() = d;
440  }
441 
442  ++derives;
443  ++points;
444  ++displacements;
445  }
446 
447 }
448 
449 /* Copy the content of m_Location into the Output. */
450 template <typename TInputMesh, typename TOutputMesh>
451 void
452 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
453 ::ComputeOutput()
454 {
455 
456  int i;
457  typename TriCell::CellAutoPointer insertCell;
458  unsigned long tripoints[3];
459  const unsigned long *tp;
460  double x;
461 
462  OutputMeshType * output = this->GetOutput();
463 
464  typedef typename OutputMeshType::CellsContainer CellsContainer;
465 
466  output->SetCells(CellsContainer::New());
467 
468  output->GetCells()->Reserve( m_NumberOfCells );
469 
470  output->SetCellsAllocationMethod( OutputMeshType::CellsAllocatedDynamicallyCellByCell );
471 
472  OutputPointsContainerPointer myPoints = output->GetPoints();
473  myPoints->Reserve(m_NumberOfNodes);
474  OutputPointsContainerIterator points = myPoints->Begin();
475 
476  InputPointsContainerPointer myLocations = m_Locations->GetPoints();
477  InputPointsContainerIterator locations = myLocations->Begin();
478 
479  InputMeshConstPointer inputMesh = this->GetInput(0);
480  InputCellsContainerConstPointer myCells = inputMesh->GetCells();
481  InputCellsContainerConstIterator cells = myCells->Begin();
482 
483  InputCellDataContainerConstPointer myCellData = inputMesh->GetCellData();
484  InputCellDataContainerConstIterator celldata = myCellData->Begin();
485 
486  i = 0;
487  for (; i<m_NumberOfNodes; i++)
488  {
489  points.Value() = locations.Value();
490  ++locations;
491  ++points;
492  }
493 
494  for (i=0; i<m_NumberOfCells; i++)
495  {
496  tp = cells.Value()->GetPointIds();
497  tripoints[0] = tp[0];
498  tripoints[1] = tp[1];
499  tripoints[2] = tp[2];
500  insertCell.TakeOwnership( new TriCell );
501  insertCell->SetPointIds(tripoints);
502  output->SetCell(i, insertCell );
503  x = celldata.Value();
504  output->SetCellData(i, (PixelType)x);
505  ++cells;
506  ++celldata;
507  }
508 
509 }
510 
511 /* Generate Data */
512 template <typename TInputMesh, typename TOutputMesh>
513 void
514 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
515 ::GenerateData()
516 {
517  this->Initialize();
518  this->SetMeshStiffness();
519 
520  while (m_Step < m_StepThreshold)
521  {
522 
523  const float progress =
524  static_cast<float>( m_Step ) /
525  static_cast<float>( m_StepThreshold );
526 
527  this->UpdateProgress( progress );
528  this->ComputeNormals();
529  this->GradientFit();
530 
531  if ( m_PotentialOn )
532  {
533  this->PotentialFit();
534  }
535 
536  this->ComputeDt();
537  this->Advance();
538  m_Step++;
539  }
540 
541  this->ComputeOutput();
542 }
543 
545 template <typename TInputMesh, typename TOutputMesh>
546 void
547 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
548 ::PotentialFit()
549 {
550 
551  PixelType max, extends[3], t, xs, ys, zs;
552  typename TInputMesh::PointType vec_for, vec_nor, vec_p, vec_1, vec_2;
553  int i, label;
554  ImageIndexType coord = {{0, 0, 0}};
555  ImageIndexType extend = {{0, 0, 0}};
556  int flag=0;
557 
558  InputPointsContainerPointer Points = m_Locations->GetPoints();
559  InputPointsContainerIterator points = Points->Begin();
560 
561  InputPointsContainerPointer myForces = m_Forces->GetPoints();
562  InputPointsContainerIterator forces = myForces->Begin();
563 
564  InputPointsContainerPointer myNormals = m_Normals->GetPoints();
565  InputPointsContainerIterator normals = myNormals->Begin();
566 
567  i = 0;
568 
569 
570  while( i < m_NumberOfNodes )
571  {
572  xs = ys = zs = 1.0;
573  vec_p = points.Value();
574 
575  coord[0] = (int) vec_p[0];
576  coord[1] = (int) vec_p[1];
577  coord[2] = (int) vec_p[2];
578 
579  if ( m_Potential->GetPixel(coord) != m_ObjectLabel )
580  {
581  xs = ys = zs = -1.0;
582  flag = 1;
583  }
584  extends[0] = vec_p[0];
585  extends[1] = vec_p[1];
586  extends[2] = vec_p[2];
587  extend[0] = (int) vec_p[0];
588  extend[1] = (int) vec_p[1];
589  extend[2] = (int) vec_p[2];
590 
591  vec_nor = normals.Value();
592 
593  max = vcl_abs(vec_nor[0]);
594 
595  //---------------------------------------------------------------------
596  // all the movement in z direction is now disabled for further test
597  //---------------------------------------------------------------------
598  if ( vcl_abs(vec_nor[1]) > max ) max = vcl_abs(vec_nor[1]);
599  if ( vcl_abs(vec_nor[2]) > max ) max = vcl_abs(vec_nor[2]);
600  if ( flag )
601  {
602  vec_1[0] = -1*vec_nor[0]/max;
603  vec_1[1] = -1*vec_nor[1]/max;
604  vec_1[2] = -1*vec_nor[2]/max;
605  }
606  else
607  {
608  vec_1[0] = vec_nor[0]/max;
609  vec_1[1] = vec_nor[1]/max;
610  vec_1[2] = vec_nor[2]/max;
611  }
612 
613  t = 0.0;
614 
615  while (t < 5.0)
616  {
617  extends[0] += vec_1[0];
618  extends[1] += vec_1[1];
619  extends[2] += vec_1[2];
620  extend[0] = (int) (extends[0]+1);
621  extend[1] = (int) (extends[1]+1);
622  extend[2] = (int) (extends[2]+1);
623  if ((extend[0] <= 0) || (extend[1] <= 0) || (extend[2] <= 0)) break;
624 
625  extend[0] = (int) (extends[0]);
626  extend[1] = (int) (extends[1]);
627  extend[2] = (int) (extends[2]);
628  if ((extend[0] >= m_ImageWidth) || (extend[1] >= m_ImageHeight) ||
629  (extend[2] >= m_ImageDepth)) break;
630 
631  label = m_Potential->GetPixel(extend);
632  if ( !flag )
633  {
634  if ( label != m_ObjectLabel ) break;
635  }
636  else if ( label == m_ObjectLabel ) break;
637 
638  t += 1.0;
639  }
640 
641  vec_2[0] = t*m_PotentialMagnitude*vec_nor[0]*xs;
642  vec_2[1] = t*m_PotentialMagnitude*vec_nor[1]*ys;
643  vec_2[2] = t*m_PotentialMagnitude*vec_nor[2]*zs;
644 
645  vec_for = forces.Value();
646  vec_for[0] += vec_2[0];
647  vec_for[1] += vec_2[1];
648  vec_for[2] += vec_2[2];
649  forces.Value() = vec_for;
650 
651  ++forces;
652  ++points;
653  ++normals;
654  ++i;
655  }
656 
657 }
658 
660 template <typename TInputMesh, typename TOutputMesh>
661 void
662 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
663 ::GradientFit()
664 {
665  typedef typename ImageIndexType::IndexValueType IndexValueType;
666 
667  ImageIndexType coord, coord2, tmp_co_1, tmp_co_2, tmp_co_3;
668  InputPointType v1, v2;
669  PixelType mag;
670 
671  typename TInputMesh::PointType vec_nor, vec_loc, vec_for, tmp_vec_1, tmp_vec_2, tmp_vec_3;
672 
673  InputPointsContainerPointer myLocations = m_Locations->GetPoints();
674  InputPointsContainerIterator locations = myLocations->Begin();
675 
676  InputPointsContainerPointer myForces = m_Forces->GetPoints();
677  InputPointsContainerIterator forces = myForces->Begin();
678 
679  InputPointDataContainerPointer myForceData = m_Forces->GetPointData();
680  InputPointDataContainerIterator forcedata = myForceData->Begin();
681 
682  InputPointsContainerPointer myNormals = m_Normals->GetPoints();
683  InputPointsContainerIterator normals = myNormals->Begin();
684 
685  /* New gradient fit method testing. */
686  while( forces != myForces->End() )
687  {
688  vec_loc = locations.Value();
689  vec_nor = normals.Value();
690 
691  coord[0] = static_cast<IndexValueType>(vec_loc[0]);
692  coord[1] = static_cast<IndexValueType>(vec_loc[1]);
693  coord[2] = static_cast<IndexValueType>(vec_loc[2]);
694 
695  coord2[0] = static_cast<IndexValueType>( vcl_ceil(vec_loc[0]) );
696  coord2[1] = static_cast<IndexValueType>( vcl_ceil(vec_loc[1]) );
697  coord2[2] = static_cast<IndexValueType>( vcl_ceil(vec_loc[2]) );
698 
699  tmp_co_1[0] = coord2[0];
700  tmp_co_1[1] = coord[1];
701  tmp_co_1[2] = coord[2];
702 
703  tmp_co_2[0] = coord[0];
704  tmp_co_2[1] = coord2[1];
705  tmp_co_2[2] = coord[2];
706 
707  tmp_co_3[0] = coord[0];
708  tmp_co_3[1] = coord[1];
709  tmp_co_3[2] = coord2[2];
710 
711  if ( (coord[0] >= 0) && (coord[1] >= 0) && (coord[2] >= 0) &&
712  (coord2[0] < m_ImageWidth) && (coord2[1] < m_ImageHeight) && (coord2[2] < m_ImageDepth) )
713  {
714  vec_for[0] = m_Gradient->GetPixel(coord)[0];
715  vec_for[1] = m_Gradient->GetPixel(coord)[1];
716  vec_for[2] = m_Gradient->GetPixel(coord)[2];
717 
718  tmp_vec_1[0] = m_Gradient->GetPixel(tmp_co_1)[0] - m_Gradient->GetPixel(coord)[0];
719  tmp_vec_1[1] = m_Gradient->GetPixel(tmp_co_1)[1] - m_Gradient->GetPixel(coord)[1];
720  tmp_vec_1[2] = m_Gradient->GetPixel(tmp_co_1)[2] - m_Gradient->GetPixel(coord)[2];
721  tmp_vec_2[0] = m_Gradient->GetPixel(tmp_co_2)[0] - m_Gradient->GetPixel(coord)[0];
722  tmp_vec_2[1] = m_Gradient->GetPixel(tmp_co_2)[1] - m_Gradient->GetPixel(coord)[1];
723  tmp_vec_2[2] = m_Gradient->GetPixel(tmp_co_2)[2] - m_Gradient->GetPixel(coord)[2];
724  tmp_vec_3[0] = m_Gradient->GetPixel(tmp_co_3)[0] - m_Gradient->GetPixel(coord)[0];
725  tmp_vec_3[1] = m_Gradient->GetPixel(tmp_co_3)[1] - m_Gradient->GetPixel(coord)[1];
726  tmp_vec_3[2] = m_Gradient->GetPixel(tmp_co_3)[2] - m_Gradient->GetPixel(coord)[2];
727 
728  vec_for[0] = vec_for[0] + (vec_loc[0]-coord[0])*tmp_vec_1[0]
729  + (vec_loc[1]-coord[1])*tmp_vec_2[0] + (vec_loc[2]-coord[2])*tmp_vec_3[0];
730  vec_for[1] = vec_for[1] + (vec_loc[1]-coord[1])*tmp_vec_2[1]
731  + (vec_loc[0]-coord[0])*tmp_vec_1[1] + (vec_loc[2]-coord[2])*tmp_vec_3[1];
732  vec_for[2] = vec_for[2] + (vec_loc[2]-coord[2])*tmp_vec_3[2]
733  + (vec_loc[1]-coord[1])*tmp_vec_2[2] + (vec_loc[0]-coord[0])*tmp_vec_1[2];
734  }
735  else
736  {
737  vec_for[0] = 0;
738  vec_for[1] = 0;
739  vec_for[2] = 0;
740  }
741 
742  mag = vec_for[0]*vec_nor[0] + vec_for[1]*vec_nor[1]+ vec_for[2]*vec_nor[2];
743 
744  vec_for[0] = m_GradientMagnitude*mag*vec_nor[0]/*num_for*/;
745  vec_for[1] = m_GradientMagnitude*mag*vec_nor[1]/*num_for*/;
746  vec_for[2] = m_GradientMagnitude*mag*vec_nor[2]/*num_for*/;
747 
748  mag = vcl_sqrt(vec_for[0]*vec_for[0] + vec_for[1]*vec_for[1]+ vec_for[2]*vec_for[2]);
749  if (mag > 0.5)
750  {
751  for (int i=0; i<3; i++)
752  {
753  vec_for[i] = (0.5 * vec_for[i])/mag;
754  }
755  }
756  forces.Value() = vec_for;
757 
758  ++forces;
759  ++forcedata;
760  ++locations;
761  ++normals;
762  }
763 }
764 
765 /* Compute normals. */
766 template <typename TInputMesh, typename TOutputMesh>
767 void
768 DeformableMesh3DFilter<TInputMesh, TOutputMesh>
769 ::ComputeNormals()
770 {
771  const unsigned long *tp;
772  InputPointType v1, v2, v3, v4, d;
773  v1.Fill(0.);
774  v2.Fill(0.);
775  v3.Fill(0.);
776  d.Fill(0.);
777 
778  double coa, cob, coc;
779  double absvec;
780 
781  InputMeshConstPointer inputMesh = this->GetInput(0);
782  InputCellsContainerConstPointer myCells = inputMesh->GetCells();
783  InputCellsContainerConstIterator cells = myCells->Begin();
784 
785  InputPointsContainerPointer myNormals = m_Normals->GetPoints();
786  InputPointsContainerIterator normals = myNormals->Begin();
787 
788  while( normals != myNormals->End() )
789  {
790  normals.Value() = d;
791  ++normals;
792  }
793 
794  while ( cells != myCells->End() )
795  {
796  tp = cells.Value()->GetPointIds();
797  ++cells;
798 
799  m_Locations->GetPoint (tp[0], &v1);
800  m_Locations->GetPoint (tp[1], &v2);
801  m_Locations->GetPoint (tp[2], &v3);
802 
803  coa = -(v1[1]*(v2[2]-v3[2]) +
804  v2[1]*(v3[2]-v1[2]) +
805  v3[1]*(v1[2]-v2[2]));
806  cob = -(v1[2] * (v2[0]-v3[0]) +
807  v2[2]*(v3[0]-v1[0]) +
808  v3[2]*(v1[0]-v2[0]));
809  coc = -(v1[0] * (v2[1]-v3[1]) +
810  v2[0]*(v3[1]-v1[1]) +
811  v3[0]*(v1[1]-v2[1]));
812 
813  absvec = -vcl_sqrt ((double) ((coa*coa) + (cob*cob) + (coc*coc)));
814 
815  assert (absvec != 0);
816 
817  v4[0] = coa/absvec;
818  v4[1] = cob/absvec;
819  v4[2] = coc/absvec;
820  m_Normals->GetPoint (tp[0], &v1);
821  m_Normals->GetPoint (tp[1], &v2);
822  m_Normals->GetPoint (tp[2], &v3);
823 
824  v1[0] += v4[0];
825  v1[1] += v4[1];
826  v1[2] += v4[2];
827 
828  v2[0] += v4[0];
829  v2[1] += v4[1];
830  v2[2] += v4[2];
831 
832  v3[0] += v4[0];
833  v3[1] += v4[1];
834  v3[2] += v4[2];
835 
836  m_Normals->SetPoint (tp[0], v1);
837  m_Normals->SetPoint (tp[1], v2);
838  m_Normals->SetPoint (tp[2], v3);
839 
840  }
841 
842  normals = myNormals->Begin();
843  while( normals != myNormals->End() )
844  {
845  v1 = normals.Value();
846 
847  absvec = vcl_sqrt((double) ((v1[0]*v1[0]) + (v1[1]*v1[1]) +
848  (v1[2]*v1[2])));
849  v1[0] = v1[0]/absvec;
850  v1[1] = v1[1]/absvec;
851  v1[2] = v1[2]/absvec;
852 
853  normals.Value() = v1;
854  ++normals;
855  }
856 }
857 
858 } /* end namespace itk. */
859 
860 #endif
Generated at Sat Feb 2 2013 23:34:37 for Orfeo Toolbox with doxygen 1.8.1.1