KWStyle - itkTetrahedronCell.txx
 
Matrix View
Description
1 /*=========================================================================
2
3   Program:   Insight Segmentation & Registration Toolkit
4   Module:    $RCSfile: itkTetrahedronCell.txx.html,v $
5   Language:  C++
6   Date:      $Date: 2006/01/17 19:15:48 $
7   Version:   $Revision: 1.4 $
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 DEF #ifndef _itkTetrahedronCell_txx
18 DEF #define _itkTetrahedronCell_txx
19 #include "itkTetrahedronCell.h"
20 #include "itkTriangleCell.h"
21 #include "vnl/algo/vnl_determinant.h"
22
23 namespace itk
24 {
25
26 /**
27  * Standard CellInterface:
28  */
29 template <typename TCellInterface>
30 void
31 TetrahedronCell< TCellInterface >
32 ::MakeCopy(CellAutoPointer & cellPointer) const
33 {
34   cellPointer.TakeOwnership( new Self );
35   cellPointer->SetPointIds(this->GetPointIds());
36 }
37
38
39 /**
40  * Standard CellInterface:
41  * Get the topological dimension of this cell.
42  */
43 template <typename TCellInterface>
44 unsigned int
45 TetrahedronCell< TCellInterface >
46 ::GetDimension(void) const
47 {
48   return Self::CellDimension;
49 }
50
51
52 /**
53  * Standard CellInterface:
54  * Get the number of points required to define the cell.
55  */
56 template <typename TCellInterface>
57 unsigned int
58 TetrahedronCell< TCellInterface >
59 ::GetNumberOfPoints(void) const
60 {
61   return Self::NumberOfPoints;
62 }  
63
64
65 /**
66  * Standard CellInterface:
67  * Get the number of boundary features of the given dimension.
68  */
69 template <typename TCellInterface>
70 typename TetrahedronCell< TCellInterface >::CellFeatureCount
71 TetrahedronCell< TCellInterface >
72 ::GetNumberOfBoundaryFeatures(int dimension) const
73 {
74   switch (dimension)
75     {
76     case 0: return GetNumberOfVertices();
77     case 1: return GetNumberOfEdges();
78     case 2: return GetNumberOfFaces();
79     default: return 0;
80     }
81 }
82
83
84 template <typename TCellInterface>
85 bool
86 TetrahedronCell< TCellInterface >
87 ::EvaluatePosition(CoordRepType* x,
88                    PointsContainer* points,
89                    CoordRepType* closestPoint,
90                    CoordRepType pcoord[3],
91                    double* minDist2,
92                    InterpolationWeightType* weights)
93 {
94   unsigned int i;
95   double rhs[PointDimension];
96   double c1[PointDimension]; 
97   double c2[PointDimension];
98   double c3[PointDimension];
99   double det;
100   double p4;
101
102   CoordRepType pcoords[3];
103   pcoords[0] = pcoords[1] = pcoords[2] = 0.0;
104
105   if(!points)
106     {
107     return false;
108     }
109
110   PointType pt1 = points->GetElement( m_PointIds[0] );
111   PointType pt2 = points->GetElement( m_PointIds[1] );
112   PointType pt3 = points->GetElement( m_PointIds[2] );
113   PointType pt4 = points->GetElement( m_PointIds[3] );
114
115   for (i=0; i<PointDimension; i++)
116     {
117     rhs[i] = x[i] - pt4[i];
118     c1[i] = pt1[i] - pt4[i];
119     c2[i] = pt2[i] - pt4[i];
120     c3[i] = pt3[i] - pt4[i];
121     }
122
123   // Create a vnl_matrix so that the determinant can be computed 
124   // for any PointDimension
125   vnl_matrix_fixed<CoordRepType,3,PointDimension> mat;
126   for(i=0;i<PointDimension;i++)
127     {
128     mat.put(0,i,c1[i]);
129     mat.put(1,i,c2[i]);
130     mat.put(2,i,c3[i]);
131     }
132
133   if ( (det = vnl_determinant(mat)) == 0.0 )
134     {
135     return false;
136     }
137
138   
139   for(i=0;i<PointDimension;i++)
140     {
141     mat.put(0,i,rhs[i]);
142     mat.put(1,i,c2[i]);
143     mat.put(2,i,c3[i]);
144     }
145
146   pcoords[0] = vnl_determinant(mat) / det;
147
148   for(i=0;i<PointDimension;i++)
149     {
150     mat.put(0,i,c1[i]);
151     mat.put(1,i,rhs[i]);
152     mat.put(2,i,c3[i]);
153     }
154
155   pcoords[1] = vnl_determinant(mat) / det;
156   
157   for(i=0;i<PointDimension;i++)
158     {
159     mat.put(0,i,c1[i]);
160     mat.put(1,i,c2[i]);
161     mat.put(2,i,rhs[i]);
162     }
163
164   pcoords[2] = vnl_determinant(mat) / det;
165
166   p4 = 1.0 - pcoords[0] - pcoords[1] - pcoords[2];
167
168   if(weights)
169     {
170     weights[0] = p4;
171     weights[1] = pcoords[0];
172     weights[2] = pcoords[1];
173     weights[3] = pcoords[2];
174     }
175
176   if(pcoord)
177     {
178     pcoord[0] = pcoords[0]; 
179     pcoord[1] = pcoords[1];
180     pcoord[2] = pcoords[2];
181     }
182
183   if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
184   pcoords[1] >= -0.001 && pcoords[1] <= 1.001 &&
185   pcoords[2] >= -0.001 && pcoords[2] <= 1.001 && p4 >= -0.001 && p4 <= 1.001 )
186     {
187     if (closestPoint)
188       {
189       for(unsigned int i=0;i <PointDimension;i++)
190         {
191         closestPoint[i] = x[i]; 
192         }
193       if(minDist2)
194         {
195         *minDist2 = 0.0; //inside tetra
196         }
197       }
198     return true; 
199     }
200   else
201     { //could easily be sped up using parametric localization - next release
202     double dist2;
203     CoordRepType  closest[PointDimension], pc[3];
204
205     if (closestPoint)
206       {
207       FaceAutoPointer triangle;
208       *minDist2 = NumericTraits< double >::max();
209       for (i=0; i<4; i++)
210         {
211         this->GetFace (i,triangle);
212         triangle->EvaluatePosition(x,points,closest,pc,&dist2,NULL);
213         
214         if ( dist2 < *minDist2 )
215           {
216           closestPoint[0] = closest[0]; 
217           closestPoint[1] = closest[1]; 
218           closestPoint[2] = closest[2];
219           *minDist2 = dist2;
220           }
221         }
222       }
223     //Just fall through to default return false;
224     }
225     return false;
226 }
227
228
229 /**
230  * Standard CellInterface:
231  * Get the boundary feature of the given dimension specified by the given
232  * cell feature Id.
233  * The Id can range from 0 to GetNumberOfBoundaryFeatures(dimension)-1.
234  */
235 template <typename TCellInterface>
236 bool
237 TetrahedronCell< TCellInterface >
238 ::GetBoundaryFeature(int dimension, CellFeatureIdentifier featureId, 
239                      CellAutoPointer & cellPointer )
240 {
241   switch (dimension)
242     {
243     case 0: 
244 IND ****{
245 IND ****VertexAutoPointer vertexPointer;
246 IND ****if( this->GetVertex(featureId,vertexPointer) )
247 IND ******{
248 IND ******TransferAutoPointer(cellPointer,vertexPointer);
249 IND ******return true;
250 IND ******}
251 IND ****else
252 IND ******{
253 IND ******cellPointer.Reset();
254 IND ******return false;
255 IND ******}
256 IND ****break;
257 IND ****}
258     case 1: 
259 IND ****{
260 IND ****EdgeAutoPointer edgePointer;
261 IND ****if( this->GetEdge(featureId,edgePointer) )
262 IND ******{
263 IND ******TransferAutoPointer(cellPointer,edgePointer);
264 IND ******return true;
265 IND ******}
266 IND ****else
267 IND ******{
268 IND ******cellPointer.Reset();
269 IND ******return false;
270 IND ******}
271 IND ****break;
272 IND ****}
273     case 2: 
274 IND ****{
275 IND ****FaceAutoPointer facePointer;
276 IND ****if( this->GetFace(featureId,facePointer) )
277 IND ******{
278 IND ******TransferAutoPointer(cellPointer,facePointer);
279 IND ******return true;
280 IND ******}
281 IND ****else
282 IND ******{
283 IND ******cellPointer.Reset();
284 IND ******return false;
285 IND ******}
286 IND ****break;
287 IND ****}
288     default: 
289 IND ****{
290 IND ****cellPointer.Reset();
291 IND ****return false;
292 IND ****}
293     }
294   return false;
295 }
296
297
298 /**
299  * Standard CellInterface:
300  * Set the point id list used by the cell.  It is assumed that the given
301  * iterator can be incremented and safely de-referenced enough times to 
302  * get all the point ids needed by the cell.
303  */
304 template <typename TCellInterface>
305 void
306 TetrahedronCell< TCellInterface >
307 ::SetPointIds(PointIdConstIterator first)
308 {
309   PointIdConstIterator ii(first);
310 SEM   for(unsigned int i=0; i < Self::NumberOfPoints ; ++i)
311     {
312     m_PointIds[i] = *ii++;
313     }
314 }
315
316
317 /**
318  * Standard CellInterface:
319  * Set the point id list used by the cell.  It is assumed that the range
320  * of iterators [first, last) contains the correct number of points needed to
321  * define the cell.  The position *last is NOT referenced, so it can safely
322  * be one beyond the end of an array or other container.
323  */
324 template <typename TCellInterface>
325 void
326 TetrahedronCell< TCellInterface >
327 ::SetPointIds(PointIdConstIterator first, PointIdConstIterator last)
328 {
329   int localId=0;
330   PointIdConstIterator ii(first);
331   
332   while(ii != last)
333     {
334     m_PointIds[localId++] = *ii++;
335     }
336 }
337
338
339 /**
340  * Standard CellInterface:
341  * Set an individual point identifier in the cell.
342  */
343 template <typename TCellInterface>
344 void
345 TetrahedronCell< TCellInterface >
346 ::SetPointId(int localId, PointIdentifier ptId)
347 {
348   m_PointIds[localId] = ptId;
349 }
350
351
352 /**
353  * Standard CellInterface:
354  * Get a begin iterator to the list of point identifiers used by the cell.
355  */
356 template <typename TCellInterface>
357 typename TetrahedronCell< TCellInterface >::PointIdIterator
358 TetrahedronCell< TCellInterface >
359 ::PointIdsBegin(void)
360 {
361   return &m_PointIds[0];
362 }
363
364
365 /**
366  * Standard CellInterface:
367  * Get a const begin iterator to the list of point identifiers used
368  * by the cell.
369  */
370 template <typename TCellInterface>
371 typename TetrahedronCell< TCellInterface >::PointIdConstIterator
372 TetrahedronCell< TCellInterface >
373 ::PointIdsBegin(void) const
374 {
375   return &m_PointIds[0];
376 }
377
378
379 /**
380  * Standard CellInterface:
381  * Get an end iterator to the list of point identifiers used by the cell.
382  */
383 template <typename TCellInterface>
384 typename TetrahedronCell< TCellInterface >::PointIdIterator
385 TetrahedronCell< TCellInterface >
386 ::PointIdsEnd(void)
387 {
388   return &m_PointIds[Self::NumberOfPoints];
389 }
390
391
392 /**
393  * Standard CellInterface:
394  * Get a const end iterator to the list of point identifiers used
395  * by the cell.
396  */
397 template <typename TCellInterface>
398 typename TetrahedronCell< TCellInterface >::PointIdConstIterator
399 TetrahedronCell< TCellInterface >
400 ::PointIdsEnd(void) const
401 {
402   return &m_PointIds[Self::NumberOfPoints];
403 }
404
405
406 /**
407  * Tetrahedron-specific:
408  * Get the number of vertices defining the tetrahedron.
409  */
410 template <typename TCellInterface>
411 typename TetrahedronCell< TCellInterface >::CellFeatureCount
412 TetrahedronCell< TCellInterface >
413 ::GetNumberOfVertices(void) const
414 {
415   return Self::NumberOfVertices;
416 }
417
418
419 /**
420  * Tetrahedron-specific:
421  * Get the number of edges defined for the tetrahedron.
422  */
423 template <typename TCellInterface>
424 typename TetrahedronCell< TCellInterface >::CellFeatureCount
425 TetrahedronCell< TCellInterface >
426 ::GetNumberOfEdges(void) const
427 {
428   return Self::NumberOfEdges;
429 }
430
431
432 /**
433  * Tetrahedron-specific:
434  * Get the number of faces defined for the tetrahedron.
435  */
436 template <typename TCellInterface>
437 typename TetrahedronCell< TCellInterface >::CellFeatureCount
438 TetrahedronCell< TCellInterface >
439 ::GetNumberOfFaces(void) const
440 {
441   return Self::NumberOfFaces;
442 }
443
444
445 /**
446  * Tetrahedron-specific:
447  * Get the vertex specified by the given cell feature Id.
448  * The Id can range from 0 to GetNumberOfVertices()-1.
449  */
450 template <typename TCellInterface>
451 bool
452 TetrahedronCell< TCellInterface >
453 ::GetVertex(CellFeatureIdentifier vertexId,VertexAutoPointer & vertexPointer )
454 {
455   VertexType * vert = new VertexType;
456   vert->SetPointId(0, m_PointIds[vertexId]);
457   vertexPointer.TakeOwnership( vert ); 
458   return true;
459 }
460
461
462 /**
463  * Tetrahedron-specific:
464  * Get the edge specified by the given cell feature Id.
465  * The Id can range from 0 to GetNumberOfEdges()-1.
466  */
467 template <typename TCellInterface>
468 bool
469 TetrahedronCell< TCellInterface >
470 ::GetEdge(CellFeatureIdentifier edgeId, EdgeAutoPointer & edgePointer )
471 {
472   EdgeType * edge = new EdgeType;
473   for(int i=0; i < EdgeType::NumberOfPoints; ++i)
474     {
475     edge->SetPointId(i, m_PointIds[ m_Edges[edgeId][i] ]);
476     }
477   edgePointer.TakeOwnership( edge ); 
478   return true;
479 }
480
481
482 EML
483 /**
484  * Tetrahedron-specific:
485  * Get the face specified by the given cell feature Id.
486  * The Id can range from 0 to GetNumberOfFaces()-1.
487  */
488 template <typename TCellInterface>
489 bool
490 TetrahedronCell< TCellInterface >
491 ::GetFace(CellFeatureIdentifier faceId, FaceAutoPointer & facePointer )
492 {
493   FaceType * face = new FaceType;
494   for(unsigned int i=0; i < FaceType::NumberOfPoints; ++i)
495     {
496     face->SetPointId(i, m_PointIds[ m_Faces[faceId][i] ]);
497     }
498   facePointer.TakeOwnership( face ); 
499   return true;
500 }
501
502
503 } // end namespace itk
504
505 #endif
506
Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:49PM
© Kitware Inc.