KWStyle - itkHexahedronCell.txx
 
Matrix View
Description
1 /*=========================================================================
2
3   Program:   Insight Segmentation & Registration Toolkit
4   Module:    $RCSfile: itkHexahedronCell.txx.html,v $
5   Language:  C++
6   Date:      $Date: 2006/01/17 19:15:36 $
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 _itkHexahedronCell_txx
18 DEF #define _itkHexahedronCell_txx
19 #include "itkHexahedronCell.h"
20 #include "vnl/vnl_matrix_fixed.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 HexahedronCell< TCellInterface >
32 ::MakeCopy(CellAutoPointer & cellPointer) const
33 {
34   cellPointer.TakeOwnership( new Self );
35   cellPointer->SetPointIds(this->GetPointIds());
36 }
37
38 /**
39  * Standard CellInterface:
40  * Get the topological dimension of this cell.
41  */
42 template <typename TCellInterface>
43 unsigned int
44 HexahedronCell< TCellInterface >
45 ::GetDimension(void) const
46 {
47   return Self::CellDimension;
48 }
49
50 /**
51  * Standard CellInterface:
52  * Get the number of points required to define the cell.
53  */
54 template <typename TCellInterface>
55 unsigned int
56 HexahedronCell< TCellInterface >
57 ::GetNumberOfPoints(void) const
58 {
59   return Self::NumberOfPoints;
60 }  
61
62 /**
63  * Standard CellInterface:
64  * Get the number of boundary features of the given dimension.
65  */
66 template <typename TCellInterface>
67 typename HexahedronCell< TCellInterface >::CellFeatureCount
68 HexahedronCell< TCellInterface >
69 ::GetNumberOfBoundaryFeatures(int dimension) const
70 {
71   switch (dimension)
72     {
73     case 0: return GetNumberOfVertices();
74     case 1: return GetNumberOfEdges();
75     case 2: return GetNumberOfFaces();
76     default: return 0;
77     }
78 }
79
80 /**
81  * Standard CellInterface:
82  * Get the boundary feature of the given dimension specified by the given
83  * cell feature Id.
84  * The Id can range from 0 to GetNumberOfBoundaryFeatures(dimension)-1.
85  */
86 template <typename TCellInterface>
87 bool
88 HexahedronCell< TCellInterface >
89 ::GetBoundaryFeature(int dimension, CellFeatureIdentifier featureId, 
90                      CellAutoPointer & cellPointer )
91 {
92   switch (dimension)
93     {
94     case 0: 
95 IND ****{
96 IND ****VertexAutoPointer vertexPointer;
97 IND ****if( this->GetVertex(featureId,vertexPointer) )
98 IND ******{
99 IND ******TransferAutoPointer(cellPointer,vertexPointer);
100 IND ******return true;
101 IND ******}
102 IND ****else
103 IND ******{
104 IND ******cellPointer.Reset();
105 IND ******return false;
106 IND ******}
107 IND ****break;
108 IND ****}
109     case 1: 
110 IND ****{
111 IND ****EdgeAutoPointer edgePointer;
112 IND ****if( this->GetEdge(featureId,edgePointer) )
113 IND ******{
114 IND ******TransferAutoPointer(cellPointer,edgePointer);
115 IND ******return true;
116 IND ******}
117 IND ****else
118 IND ******{
119 IND ******cellPointer.Reset();
120 IND ******return false;
121 IND ******}
122 IND ****break;
123 IND ****}
124     case 2: 
125 IND ****{
126 IND ****FaceAutoPointer facePointer;
127 IND ****if( this->GetFace(featureId,facePointer) )
128 IND ******{
129 IND ******TransferAutoPointer(cellPointer,facePointer);
130 IND ******return true;
131 IND ******}
132 IND ****else
133 IND ******{
134 IND ******cellPointer.Reset();
135 IND ******return false;
136 IND ******}
137 IND ****break;
138 IND ****}
139     default: 
140 IND ****{
141 IND ****cellPointer.Reset();
142 IND ****return false;
143 IND ****}
144     }
145   return false;
146 }
147
148 /**
149  * Standard CellInterface:
150  * Set the point id list used by the cell.  It is assumed that the given
151  * iterator can be incremented and safely de-referenced enough times to 
152  * get all the point ids needed by the cell.
153  */
154 template <typename TCellInterface>
155 void
156 HexahedronCell< TCellInterface >
157 ::SetPointIds(PointIdConstIterator first)
158 {
159   PointIdConstIterator ii(first);
160 SEM   for(int i=0; i < Self::NumberOfPoints ; ++i)
161     {   
162     m_PointIds[i] = *ii++;
163     }
164 }
165
166 /**
167  * Standard CellInterface:
168  * Set the point id list used by the cell.  It is assumed that the range
169  * of iterators [first, last) contains the correct number of points needed to
170  * define the cell.  The position *last is NOT referenced, so it can safely
171  * be one beyond the end of an array or other container.
172  */
173 template <typename TCellInterface>
174 void
175 HexahedronCell< TCellInterface >
176 ::SetPointIds(PointIdConstIterator first, PointIdConstIterator last)
177 {
178   int localId=0;
179   PointIdConstIterator ii(first);
180   
181   while(ii != last)
182     {
183     m_PointIds[localId++] = *ii++;
184     }
185 }
186
187 /**
188  * Standard CellInterface:
189  * Set an individual point identifier in the cell.
190  */
191 template <typename TCellInterface>
192 void
193 HexahedronCell< TCellInterface >
194 ::SetPointId(int localId, PointIdentifier ptId)
195 {
196   m_PointIds[localId] = ptId;
197 }
198
199 /**
200  * Standard CellInterface:
201  * Get a begin iterator to the list of point identifiers used by the cell.
202  */
203 template <typename TCellInterface>
204 typename HexahedronCell< TCellInterface >::PointIdIterator
205 HexahedronCell< TCellInterface >
206 ::PointIdsBegin(void)
207 {
208   return &m_PointIds[0];
209 }
210
211 /**
212  * Standard CellInterface:
213  * Get a const begin iterator to the list of point identifiers used
214  * by the cell.
215  */
216 template <typename TCellInterface>
217 typename HexahedronCell< TCellInterface >::PointIdConstIterator
218 HexahedronCell< TCellInterface >
219 ::PointIdsBegin(void) const
220 {
221   return &m_PointIds[0];
222 }
223
224 /**
225  * Standard CellInterface:
226  * Get an end iterator to the list of point identifiers used by the cell.
227  */
228 template <typename TCellInterface>
229 typename HexahedronCell< TCellInterface >::PointIdIterator
230 HexahedronCell< TCellInterface >
231 ::PointIdsEnd(void)
232 {
233   return &m_PointIds[Self::NumberOfPoints];
234 }
235
236 /**
237  * Standard CellInterface:
238  * Get a const end iterator to the list of point identifiers used
239  * by the cell.
240  */
241 template <typename TCellInterface>
242 typename HexahedronCell< TCellInterface >::PointIdConstIterator
243 HexahedronCell< TCellInterface >
244 ::PointIdsEnd(void) const
245 {
246   return &m_PointIds[Self::NumberOfPoints];
247 }
248
249 /**
250  * Hexahedron-specific:
251  * Get the number of vertices defining the hexahedron.
252  */
253 template <typename TCellInterface>
254 typename HexahedronCell< TCellInterface >::CellFeatureCount
255 HexahedronCell< TCellInterface >
256 ::GetNumberOfVertices(void) const
257 {
258   return Self::NumberOfVertices;
259 }
260
261 /**
262  * Hexahedron-specific:
263  * Get the number of edges defined for the hexahedron.
264  */
265 template <typename TCellInterface>
266 typename HexahedronCell< TCellInterface >::CellFeatureCount
267 HexahedronCell< TCellInterface >
268 ::GetNumberOfEdges(void) const
269 {
270   return Self::NumberOfEdges;
271 }
272
273 /**
274  * Hexahedron-specific:
275  * Get the number of faces defined for the hexahedron.
276  */
277 template <typename TCellInterface>
278 typename HexahedronCell< TCellInterface >::CellFeatureCount
279 HexahedronCell< TCellInterface >
280 ::GetNumberOfFaces(void) const
281 {
282   return Self::NumberOfFaces;
283 }
284
285 /**
286  * Hexahedron-specific:
287  * Get the vertex specified by the given cell feature Id.
288  * The Id can range from 0 to GetNumberOfVertices()-1.
289  */
290 template <typename TCellInterface>
291 bool
292 HexahedronCell< TCellInterface >
293 ::GetVertex(CellFeatureIdentifier vertexId,VertexAutoPointer & vertexPointer )
294 {
295   VertexType * vert = new VertexType;
296   vert->SetPointId(0, m_PointIds[vertexId]);
297   vertexPointer.TakeOwnership( vert );
298   return true;  
299 }
300
301 /**
302  * Hexahedron-specific:
303  * Get the edge specified by the given cell feature Id.
304  * The Id can range from 0 to GetNumberOfEdges()-1.
305  */
306 template <typename TCellInterface>
307 bool
308 HexahedronCell< TCellInterface >
309 ::GetEdge(CellFeatureIdentifier edgeId, EdgeAutoPointer & edgePointer )
310 {
311   EdgeType * edge = new EdgeType;
312   for(int i=0; i < EdgeType::NumberOfPoints; ++i)
313     {
314     edge->SetPointId(i, m_PointIds[ m_Edges[edgeId][i] ]);
315     }
316   edgePointer.TakeOwnership( edge ); 
317   return true;
318 }
319
320
321 /**
322  * Hexahedron-specific:
323  * Get the face specified by the given cell feature Id.
324  * The Id can range from 0 to GetNumberOfFaces()-1.
325  */
326 template <typename TCellInterface>
327 bool
328 HexahedronCell< TCellInterface >
329 ::GetFace(CellFeatureIdentifier faceId, FaceAutoPointer & facePointer )
330 {
331   FaceType * face = new FaceType;
332   for(unsigned int i=0; i < FaceType::NumberOfPoints; ++i)
333     {
334     face->SetPointId(i, m_PointIds[ m_Faces[faceId][i] ]);
335     }
336   facePointer.TakeOwnership( face ); 
337   return true;
338 }
339
340 /** Evaluate the position inside the cell */
341 template <typename TCellInterface>
342 bool
343 HexahedronCell< TCellInterface >
344 ::EvaluatePosition(CoordRepType* x,
345                    PointsContainer* points,
346                    CoordRepType* closestPoint,
347                    CoordRepType pcoord[3],
348                    double* dist2,
349                    InterpolationWeightType* weight)
350 {
351   static const int ITK_HEX_MAX_ITERATION=10;
352   static const double ITK_HEX_CONVERGED=1.e-03;
353   static const double ITK_DIVERGED = 1.e6;
354
355   int iteration, converged;
356   double  params[3];
357   double  fcol[3], rcol[3], scol[3], tcol[3];
358   int i, j;
359   double  d;
360   PointType pt;
361   CoordRepType derivs[24];
362   InterpolationWeightType weights[8];
363
364   //  set initial position for Newton's method
365   int subId = 0;
366   CoordRepType pcoords[3];
367   pcoords[0] = pcoords[1] = pcoords[2] = params[0] = params[1] = params[2]=0.5;
368
369   //  enter iteration loop
370   for (iteration=converged=0;
371        !converged && (iteration < ITK_HEX_MAX_ITERATION);  iteration++) 
372     {
373     //  calculate element interpolation functions and derivatives
374     this->InterpolationFunctions(pcoords, weights);
375     this->InterpolationDerivs(pcoords, derivs);
376
377     //  calculate newton functions
378     for (i=0; i<3; i++) 
379       {
380       fcol[i] = rcol[i] = scol[i] = tcol[i] = 0.0;
381       }
382     for (i=0; i<8; i++)
383       {
384       pt = points->GetElement( m_PointIds[i] );
385       for (j=0; j<3; j++)
386         {
387         fcol[j] += pt[j] * weights[i];
388         rcol[j] += pt[j] * derivs[i];
389         scol[j] += pt[j] * derivs[i+8];
390         tcol[j] += pt[j] * derivs[i+16];
391         }
392       }
393
394     for (i=0; i<3; i++)
395       {
396       fcol[i] -= x[i];
397       }
398
399     //  compute determinants and generate improvements
400     vnl_matrix_fixed<CoordRepType,3,PointDimension> mat;
401     for(i=0;i<PointDimension;i++)
402       {
403       mat.put(0,i,rcol[i]);
404       mat.put(1,i,scol[i]);
405       mat.put(2,i,tcol[i]);
406      }
407
408     d = vnl_determinant(mat);
409     //d=vtkMath::Determinant3x3(rcol,scol,tcol);
410     if ( fabs(d) < 1.e-20) 
411       {
412       return false;
413       }
414
415     vnl_matrix_fixed<CoordRepType,3,PointDimension> mat1;
416     for(i=0;i<PointDimension;i++)
417       {
418       mat1.put(0,i,fcol[i]);
419       mat1.put(1,i,scol[i]);
420       mat1.put(2,i,tcol[i]);
421      }
422
423     vnl_matrix_fixed<CoordRepType,3,PointDimension> mat2;
424     for(i=0;i<PointDimension;i++)
425       {
426       mat2.put(0,i,rcol[i]);
427       mat2.put(1,i,fcol[i]);
428       mat2.put(2,i,tcol[i]);
429      }
430
431     vnl_matrix_fixed<CoordRepType,3,PointDimension> mat3;
432     for(i=0;i<PointDimension;i++)
433       {
434       mat3.put(0,i,rcol[i]);
435       mat3.put(1,i,scol[i]);
436       mat3.put(2,i,fcol[i]);
437      }
438
439     pcoords[0] = params[0] - vnl_determinant(mat1) / d;
440     pcoords[1] = params[1] - vnl_determinant(mat2) / d;
441     pcoords[2] = params[2] - vnl_determinant(mat3) / d;
442
443     if(pcoord)
444       {
445       pcoord[0] = pcoords[0]; 
446       pcoord[1] = pcoords[1];
447       pcoord[2] = pcoords[2];
448       }
449
450     //  check for convergence
451     if ( ((fabs(pcoords[0]-params[0])) < ITK_HEX_CONVERGED) &&
452          ((fabs(pcoords[1]-params[1])) < ITK_HEX_CONVERGED) &&
453          ((fabs(pcoords[2]-params[2])) < ITK_HEX_CONVERGED) )
454       {
455       converged = 1;
456       }
457
458     // Test for bad divergence (S.Hirschberg 11.12.2001)
459     else if ((fabs(pcoords[0]) > ITK_DIVERGED) || 
460 IND *************(fabs(pcoords[1]) > ITK_DIVERGED) || 
461 IND *************(fabs(pcoords[2]) > ITK_DIVERGED))
462       {
463       return -1;
464       }
465
466     //  if not converged, repeat
467     else 
468       {
469       params[0] = pcoords[0];
470       params[1] = pcoords[1];
471       params[2] = pcoords[2];
472       }
473     }
474
475   //  if not converged, set the parametric coordinates to arbitrary values
476   //  outside of element
477   if ( !converged )
478     {
479     return false;
480     }
481
482   this->InterpolationFunctions(pcoords, weights);
483
484   if(weight)
485     {
486     for(unsigned int i=0;i<8;i++)
487       {
488       weight[i] = weights[i];
489       }
490     }
491
492   if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001 &&
493   pcoords[1] >= -0.001 && pcoords[1] <= 1.001 &&
494   pcoords[2] >= -0.001 && pcoords[2] <= 1.001 )
495     {
496     if (closestPoint)
497       {
498       closestPoint[0] = x[0]; closestPoint[1] = x[1]; closestPoint[2] = x[2];
499       *dist2 = 0.0; //inside hexahedron
500       }
501     return true;
502     }
503   else
504     {
505     CoordRepType pc[3], w[8];
506     if (closestPoint)
507       {
508       for (i=0; i<3; i++) //only approximate, not really true for warped hexa
509         {
510         if (pcoords[i] < 0.0)
511           {
512           pc[i] = 0.0;
513           }
514         else if (pcoords[i] > 1.0)
515           {
516           pc[i] = 1.0;
517           }
518         else
519           {
520           pc[i] = pcoords[i];
521           }
522         }
523 LEN       this->EvaluateLocation(subId, points, pc, closestPoint, (CoordRepType *)w);
524
525       *dist2 = 0;
526       for(unsigned int i=0;i<3;i++)
527         {
528         *dist2 += (closestPoint[i]-x[i])*(closestPoint[i]-x[i]);
529         }
530       }
531     return false;
532     }
533 }
534
535 /** Compute iso-parametric interpolation functions */
536 template <typename TCellInterface>
537 void 
538 HexahedronCell< TCellInterface >
539 ::InterpolationFunctions(CoordRepType pcoords[3], CoordRepType sf[8])
540 {
541   double rm, sm, tm;
542
543   rm = 1. - pcoords[0];
544   sm = 1. - pcoords[1];
545   tm = 1. - pcoords[2];
546
547   sf[0] = rm*sm*tm;
548   sf[1] = pcoords[0]*sm*tm;
549   sf[2] = pcoords[0]*pcoords[1]*tm;
550   sf[3] = rm*pcoords[1]*tm;
551   sf[4] = rm*sm*pcoords[2];
552   sf[5] = pcoords[0]*sm*pcoords[2];
553   sf[6] = pcoords[0]*pcoords[1]*pcoords[2];
554   sf[7] = rm*pcoords[1]*pcoords[2];
555 }
556
557 /** Compute iso-parametric interpolation functions */
558 template <typename TCellInterface>
559 void 
560 HexahedronCell< TCellInterface >
561 ::InterpolationDerivs(CoordRepType pcoords[3], CoordRepType derivs[24])
562 {
563   double rm, sm, tm;
564
565   rm = 1. - pcoords[0];
566   sm = 1. - pcoords[1];
567   tm = 1. - pcoords[2];
568
569   // r-derivatives
570   derivs[0] = -sm*tm;
571   derivs[1] = sm*tm;
572   derivs[2] = pcoords[1]*tm;
573   derivs[3] = -pcoords[1]*tm;
574   derivs[4] = -sm*pcoords[2];
575   derivs[5] = sm*pcoords[2];
576   derivs[6] = pcoords[1]*pcoords[2];
577   derivs[7] = -pcoords[1]*pcoords[2];
578
579   // s-derivatives
580   derivs[8] = -rm*tm;
581   derivs[9] = -pcoords[0]*tm;
582   derivs[10] = pcoords[0]*tm;
583   derivs[11] = rm*tm;
584   derivs[12] = -rm*pcoords[2];
585   derivs[13] = -pcoords[0]*pcoords[2];
586   derivs[14] = pcoords[0]*pcoords[2];
587   derivs[15] = rm*pcoords[2];
588
589   // t-derivatives
590   derivs[16] = -rm*sm;
591   derivs[17] = -pcoords[0]*sm;
592   derivs[18] = -pcoords[0]*pcoords[1];
593   derivs[19] = -rm*pcoords[1];
594   derivs[20] = rm*sm;
595   derivs[21] = pcoords[0]*sm;
596   derivs[22] = pcoords[0]*pcoords[1];
597   derivs[23] = rm*pcoords[1];
598 }
599
600 /** Evaluate the location inside the cell */
601 template <typename TCellInterface>
602 void 
603 HexahedronCell< TCellInterface >
604 LEN ::EvaluateLocation(int& itkNotUsed(subId),PointsContainer* points, CoordRepType pcoords[3],
605 IND *******************CoordRepType x[3], CoordRepType *weights)
606 {
607   int i, j;
608   PointType pt;
609
610   this->InterpolationFunctions(pcoords, weights);
611
612   x[0] = x[1] = x[2] = 0.0;
613   for (i=0; i<8; i++)
614     {
615     pt = points->GetElement( m_PointIds[i] );
616      
617     for (j=0; j<3; j++)
618       {
619       x[j] += pt[j] * weights[i];
620       }
621     }
622 }
623
624
625 } // end namespace itk
626
627 #endif
628
Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:12PM
© Kitware Inc.