Orfeo Toolbox  4.0
itkSimplexMeshGeometry.cxx
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #include "itkSimplexMeshGeometry.h"
19 
20 #include "vxl_version.h"
21 #if VXL_VERSION_DATE_FULL > 20040406
22 #include "vnl/vnl_cross.h"
23 #define cross_3d vnl_cross_3d
24 #endif
25 
26 namespace itk
27 {
30 {
31  double c = 1.0 / 3.0;
32  PointType p;
33 
34  p.Fill(0.0);
35 
36  pos.Fill(0);
37  oldPos.Fill(0);
38  referenceMetrics.Fill(c);
39  eps.Fill(c);
40  normal.Fill(0);
41  externalForce.Fill(0);
42  internalForce.Fill(0);
43  circleRadius = 0;
44  circleCenter.Fill(0);
45  sphereRadius = 0;
46  distance = 0;
47  phi = 0;
48 
49  neighborIndices.Fill( NumericTraits< IdentifierType >::max() );
50  neighbors.Fill(p);
51  meanCurvature = c;
52 
53  neighborSet = NULL;
54 }
55 
58 {
59  delete this->neighborSet;
60  this->neighborSet = NULL;
61 }
62 
63 void
66 {
67  VectorType b, c, cXb, tmp;
68 
69  //compute the circum circle (center and radius)
70  b = this->neighbors[2] - this->neighbors[0];
71  c = this->neighbors[1] - this->neighbors[0];
72 
73  cXb.SetVnlVector( cross_3d< double >( c.GetVnlVector(), b.GetVnlVector() ) );
74 
76  * cross_3d< double >( cXb.GetVnlVector(), c.GetVnlVector() )
77  + c.GetSquaredNorm()
78  * cross_3d< double >( b.GetVnlVector(), cXb.GetVnlVector() ) );
79 
80  double cXbSquaredNorm = 2 * cXb.GetSquaredNorm();
81 
82  circleRadius = tmp.GetNorm() / ( cXbSquaredNorm );
83  tmp[0] /= ( cXbSquaredNorm );
84  tmp[1] /= ( cXbSquaredNorm );
85  tmp[2] /= ( cXbSquaredNorm );
86  circleCenter = this->neighbors[0] + tmp;
87 
88  // Compute the circum sphere (center and radius) of a point
89  VectorType d, dXc, bXd, sphereTmp, denom;
90 
91  d = pos - this->neighbors[0];
92  dXc.SetVnlVector( cross_3d< double >( d.GetVnlVector(), c.GetVnlVector() ) );
93  bXd.SetVnlVector( cross_3d< double >( b.GetVnlVector(), d.GetVnlVector() ) );
94 
95  sphereTmp.SetVnlVector( d.GetSquaredNorm() * cXb.GetVnlVector()
96  + b.GetSquaredNorm() * dXc.GetVnlVector()
97  + c.GetSquaredNorm() * bXd.GetVnlVector() );
98 
99  double val = 2 * ( c[0] * ( b[1] * d[2] - b[2] * d[1] )
100  - c[1] * ( b[0] * d[2] - b[2] * d[0] )
101  + c[2] * ( b[0] * d[1] - b[1] * d[0] ) );
102 
103  // fix for points which lay on their neighbors plane
104  // necessary ??
105  if ( val == 0 )
106  {
107  val = 1; // itkAssertInDebugAndIgnoreInReleaseMacro (val != 0 );
108  }
109 
110  sphereRadius = sphereTmp.GetNorm() / val;
111 
112  if ( sphereRadius < 0 )
113  {
114  sphereRadius = -1 * sphereRadius;
115  }
116 }
117 
118 void
121 {
122  for( unsigned int i = 0; i < 3; i++ )
123  {
124  this->neighborIndices[i] = input.neighborIndices[i];
125  this->neighbors[i] = input.neighbors[i];
126  }
127  this->meanCurvature = input.meanCurvature;
128  this->pos = input.pos;
129  this->oldPos = input.oldPos;
130  this->eps = input.eps;
131  this->referenceMetrics = input.referenceMetrics;
132  this->normal = input.normal;
133  this->externalForce = input.externalForce;
134  this->internalForce = input.internalForce;
135  this->closestAttractor = input.closestAttractor;
137  this->circleRadius = input.circleRadius;
138  this->circleCenter = input.circleCenter;
139  this->sphereRadius = input.sphereRadius;
140  // this->sphereCenter = input.sphereCenter;
141  this->distance = input.distance;
142  this->phi = input.phi;
143  this->multiplier = input.multiplier;
144  this->forceIndex = input.forceIndex;
145  this->CopyNeigborSet( input.neighborSet );
146 }
147 
148 void
151 {
152  delete this->neighborSet;
153  if( nset )
154  {
155  this->neighborSet = new NeighborSetType;
156  this->neighborSet->insert( nset->begin(), nset->end() );
157  }
158  else
159  {
160  this->neighborSet = NULL;
161  }
162 }
163 
164 } // end namespace itk

Generated at Sat Mar 8 2014 15:33:32 for Orfeo Toolbox with doxygen 1.8.3.1