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 #ifndef itkTriangleCell_hxx
19 #define itkTriangleCell_hxx
20 #include "itkTriangleCell.h"
21 #include "vnl/algo/vnl_determinant.h"
22
23 namespace itk
24 {
25 /**
26 * Standard CellInterface:
27 */
28 template< typename TCellInterface >
29 void
30 TriangleCell< TCellInterface >
MakeCopy(CellAutoPointer & cellPointer) const31 ::MakeCopy(CellAutoPointer & cellPointer) const
32 {
33 cellPointer.TakeOwnership(new Self);
34 cellPointer->SetPointIds( this->GetPointIds() );
35 }
36
37 /**
38 * Standard CellInterface:
39 * Get the topological dimension of this cell.
40 */
41 template< typename TCellInterface >
42 unsigned int
43 TriangleCell< TCellInterface >
GetDimension() const44 ::GetDimension() const
45 {
46 return Self::CellDimension;
47 }
48
49 /**
50 * Standard CellInterface:
51 * Get the number of points required to define the cell.
52 */
53 template< typename TCellInterface >
54 unsigned int
55 TriangleCell< TCellInterface >
GetNumberOfPoints() const56 ::GetNumberOfPoints() const
57 {
58 return Self::NumberOfPoints;
59 }
60
61 /**
62 * Standard CellInterface:
63 * Get the number of boundary features of the given dimension.
64 */
65 template< typename TCellInterface >
66 typename TriangleCell< TCellInterface >::CellFeatureCount
67 TriangleCell< TCellInterface >
GetNumberOfBoundaryFeatures(int dimension) const68 ::GetNumberOfBoundaryFeatures(int dimension) const
69 {
70 switch ( dimension )
71 {
72 case 0:
73 return GetNumberOfVertices();
74 case 1:
75 return GetNumberOfEdges();
76 default:
77 return 0;
78 }
79 }
80
81 /**
82 * Standard CellInterface:
83 * Get the boundary feature of the given dimension specified by the given
84 * cell feature Id.
85 * The Id can range from 0 to GetNumberOfBoundaryFeatures(dimension)-1.
86 */
87 template< typename TCellInterface >
88 bool
89 TriangleCell< TCellInterface >
GetBoundaryFeature(int dimension,CellFeatureIdentifier featureId,CellAutoPointer & cellPointer)90 ::GetBoundaryFeature(int dimension, CellFeatureIdentifier featureId,
91 CellAutoPointer & cellPointer)
92 {
93 switch ( dimension )
94 {
95 case 0:
96 {
97 VertexAutoPointer vertexPointer;
98 if ( this->GetVertex(featureId, vertexPointer) )
99 {
100 TransferAutoPointer(cellPointer, vertexPointer);
101 return true;
102 }
103 break;
104 }
105 case 1:
106 {
107 EdgeAutoPointer edgePointer;
108 if ( this->GetEdge(featureId, edgePointer) )
109 {
110 TransferAutoPointer(cellPointer, edgePointer);
111 return true;
112 }
113 break;
114 }
115 default:
116 break; //just fall through and return false
117 }
118 cellPointer.Reset();
119 return false;
120 }
121
122 /**
123 * Standard CellInterface:
124 * Set the point id list used by the cell. It is assumed that the given
125 * iterator can be incremented and safely de-referenced enough times to
126 * get all the point ids needed by the cell.
127 */
128 template< typename TCellInterface >
129 void
130 TriangleCell< TCellInterface >
SetPointIds(PointIdConstIterator first)131 ::SetPointIds(PointIdConstIterator first)
132 {
133 PointIdConstIterator ii(first);
134
135 for ( unsigned int i = 0; i < NumberOfPoints; ++i, ++ii )
136 {
137 m_PointIds[i] = *ii;
138 }
139 }
140
141 /**
142 * Standard CellInterface:
143 * Set the point id list used by the cell. It is assumed that the range
144 * of iterators [first, last) contains the correct number of points needed to
145 * define the cell. The position *last is NOT referenced, so it can safely
146 * be one beyond the end of an array or other container.
147 */
148 template< typename TCellInterface >
149 void
150 TriangleCell< TCellInterface >
SetPointIds(PointIdConstIterator first,PointIdConstIterator last)151 ::SetPointIds(PointIdConstIterator first, PointIdConstIterator last)
152 {
153 unsigned int localId = 0;
154 PointIdConstIterator ii(first);
155
156 while ( ( ii != last ) && ( localId < NumberOfPoints ) )
157 {
158 m_PointIds[localId++] = *ii++;
159 }
160 }
161
162 /**
163 * Standard CellInterface:
164 * Set an individual point identifier in the cell.
165 */
166 template< typename TCellInterface >
167 void
168 TriangleCell< TCellInterface >
SetPointId(int localId,PointIdentifier ptId)169 ::SetPointId(int localId, PointIdentifier ptId)
170 {
171 m_PointIds[localId] = ptId;
172 }
173
174 /**
175 * Standard CellInterface:
176 * Get a begin iterator to the list of point identifiers used by the cell.
177 */
178 template< typename TCellInterface >
179 typename TriangleCell< TCellInterface >::PointIdIterator
180 TriangleCell< TCellInterface >
PointIdsBegin()181 ::PointIdsBegin()
182 {
183 return &m_PointIds[0];
184 }
185
186 /**
187 * Standard CellInterface:
188 * Get a const begin iterator to the list of point identifiers used
189 * by the cell.
190 */
191 template< typename TCellInterface >
192 typename TriangleCell< TCellInterface >::PointIdConstIterator
193 TriangleCell< TCellInterface >
PointIdsBegin() const194 ::PointIdsBegin() const
195 {
196 return &m_PointIds[0];
197 }
198
199 /**
200 * Standard CellInterface:
201 * Get an end iterator to the list of point identifiers used by the cell.
202 */
203 template< typename TCellInterface >
204 typename TriangleCell< TCellInterface >::PointIdIterator
205 TriangleCell< TCellInterface >
PointIdsEnd()206 ::PointIdsEnd()
207 {
208 return &m_PointIds[Self::NumberOfPoints - 1] + 1;
209 }
210
211 /**
212 * Standard CellInterface:
213 * Get a const end iterator to the list of point identifiers used
214 * by the cell.
215 */
216 template< typename TCellInterface >
217 typename TriangleCell< TCellInterface >::PointIdConstIterator
218 TriangleCell< TCellInterface >
PointIdsEnd() const219 ::PointIdsEnd() const
220 {
221 return &m_PointIds[Self::NumberOfPoints - 1] + 1;
222 }
223
224 /**
225 * Triangle-specific:
226 * Get the number of vertices defining the triangle.
227 */
228 template< typename TCellInterface >
229 typename TriangleCell< TCellInterface >::CellFeatureCount
230 TriangleCell< TCellInterface >
GetNumberOfVertices() const231 ::GetNumberOfVertices() const
232 {
233 return Self::NumberOfVertices;
234 }
235
236 /**
237 * Triangle-specific:
238 * Get the number of edges defined for the triangle.
239 */
240 template< typename TCellInterface >
241 typename TriangleCell< TCellInterface >::CellFeatureCount
242 TriangleCell< TCellInterface >
GetNumberOfEdges() const243 ::GetNumberOfEdges() const
244 {
245 return Self::NumberOfEdges;
246 }
247
248 /**
249 * Triangle-specific:
250 * Get the vertex specified by the given cell feature Id.
251 * The Id can range from 0 to GetNumberOfVertices()-1.
252 */
253 template< typename TCellInterface >
254 bool
255 TriangleCell< TCellInterface >
GetVertex(CellFeatureIdentifier vertexId,VertexAutoPointer & vertexPointer)256 ::GetVertex(CellFeatureIdentifier vertexId, VertexAutoPointer & vertexPointer)
257 {
258 auto * vert = new VertexType;
259
260 vert->SetPointId(0, m_PointIds[vertexId]);
261 vertexPointer.TakeOwnership(vert);
262 return true;
263 }
264
265 /**
266 * Triangle-specific:
267 * Get the edge specified by the given cell feature Id.
268 * The Id can range from 0 to GetNumberOfEdges()-1.
269 */
270 template< typename TCellInterface >
271 bool
272 TriangleCell< TCellInterface >
GetEdge(CellFeatureIdentifier edgeId,EdgeAutoPointer & edgePointer)273 ::GetEdge(CellFeatureIdentifier edgeId, EdgeAutoPointer & edgePointer)
274 {
275 auto * edge = new EdgeType;
276
277 for ( unsigned int i = 0; i < EdgeType::NumberOfPoints; ++i )
278 {
279 edge->SetPointId(i, m_PointIds[m_Edges[edgeId][i]]);
280 }
281 edgePointer.TakeOwnership(edge);
282 return true;
283 }
284
285 /** Compute distance to finite line. Returns parametric coordinate t
286 * and point location on line. */
287 template< typename TCellInterface >
288 double
289 TriangleCell< TCellInterface >
DistanceToLine(PointType x,PointType p1,PointType p2,double & t,CoordRepType * closestPoint)290 ::DistanceToLine(PointType x, PointType p1, PointType p2,
291 double & t, CoordRepType *closestPoint)
292 {
293 // convert from CoordRepType * to PointType:
294 PointType temp(closestPoint);
295 // for (unsigned int i = 0; i < PointDimension; i++)
296 // {
297 // temp[i] = closestPoint[i];
298 // }
299
300 // Compute the squared distance to the line:
301 const double distance2 = this->DistanceToLine (x, p1, p2, t, temp);
302
303 // convert from PointType to CoordRepType * :
304 for ( unsigned int j = 0; j < PointDimension; j++ )
305 {
306 closestPoint[j] = temp[j];
307 }
308
309 return distance2;
310 }
311
312 template< typename TCellInterface >
313 double
314 TriangleCell< TCellInterface >
DistanceToLine(PointType x,PointType p1,PointType p2,double & t,PointType & closestPoint)315 ::DistanceToLine(PointType x, PointType p1, PointType p2,
316 double & t, PointType & closestPoint)
317 {
318 VectorType v21 = p2 - p1;
319 //
320 // Get parametric location
321 //
322 double num(0);
323 double denom(0);
324
325 for ( unsigned int i = 0; i < PointDimension; i++ )
326 {
327 num += static_cast< double >( v21[i] * ( x[i] - p1[i] ) );
328 denom += static_cast< double >( v21[i] * v21[i] );
329 }
330
331 // trying to avoid an expensive fabs
332 double tolerance = 1.e-05 * num;
333 if ( tolerance < 0.0 )
334 {
335 tolerance = -tolerance;
336 }
337 if ( ( -tolerance < denom ) && ( denom < tolerance ) ) //numerically bad!
338 {
339 closestPoint = p1; //arbitrary, point is (numerically) far away
340 }
341 //
342 // If parametric coordinate is within 0<=p<=1, then the point is closest to
343 // the line. Otherwise, it's closest to a point at the end of the line.
344 //
345 else if ( ( t = num / denom ) < 0.0 )
346 {
347 closestPoint = p1;
348 }
349 else if ( t > 1.0 )
350 {
351 closestPoint = p2;
352 }
353 else
354 {
355 closestPoint = p1 + v21 * t;
356 }
357
358 return static_cast< double >( closestPoint.SquaredEuclideanDistanceTo(x) );
359 }
360
361 template< typename TCellInterface >
362 typename TriangleCell< TCellInterface >::CoordRepType
ComputeArea(PointsContainer * iPoints)363 TriangleCell< TCellInterface >::ComputeArea(PointsContainer *iPoints)
364 {
365 PointType p[3];
366
367 for ( unsigned int i = 0; i < NumberOfPoints; ++i )
368 {
369 p[i] = iPoints->GetElement(m_PointIds[i]);
370 }
371
372 CoordRepType a = p[1].EuclideanDistanceTo(p[2]);
373 CoordRepType b = p[0].EuclideanDistanceTo(p[2]);
374 CoordRepType c = p[1].EuclideanDistanceTo(p[0]);
375
376 CoordRepType s = 0.5 * ( a + b + c );
377 return std::sqrt( s * ( s - a ) * ( s - b ) * ( s - c ) );
378 }
379
380 template< typename TCellInterface >
381 typename TriangleCell< TCellInterface >::PointType
ComputeBarycenter(CoordRepType * iWeights,PointsContainer * iPoints)382 TriangleCell< TCellInterface >::ComputeBarycenter(
383 CoordRepType *iWeights, PointsContainer *iPoints)
384 {
385 PointType p[3];
386 CoordRepType sum_weights(0.);
387 unsigned int i(0);
388
389 for (; i < 3; i++ )
390 {
391 sum_weights += iWeights[i];
392 p[i] = iPoints->GetElement(m_PointIds[i]);
393 }
394
395 PointType oP;
396
397 if ( sum_weights != 0. )
398 {
399 oP.Fill(0.);
400 for ( i = 0; i < 3; i++ )
401 {
402 oP += p[i].GetVectorFromOrigin() * iWeights[i] / sum_weights;
403 }
404 }
405 else
406 {
407 oP = p[0];
408 }
409 return oP;
410 }
411
412 template< typename TCellInterface >
413 typename TriangleCell< TCellInterface >::PointType
ComputeCenterOfGravity(PointsContainer * iPoints)414 TriangleCell< TCellInterface >::ComputeCenterOfGravity(
415 PointsContainer *iPoints)
416 {
417 std::vector< CoordRepType > weights(3, 1. / 3.);
418 return ComputeBarycenter(& weights[0], iPoints);
419 }
420
421 template< typename TCellInterface >
422 typename TriangleCell< TCellInterface >::PointType
ComputeCircumCenter(PointsContainer * iPoints)423 TriangleCell< TCellInterface >::ComputeCircumCenter(
424 PointsContainer *iPoints)
425 {
426 std::vector< CoordRepType > weights(3, 0.);
427
428 PointType p[3];
429 unsigned int i;
430
431 for ( i = 0; i < 3; i++ )
432 {
433 p[i] = iPoints->GetElement(m_PointIds[i]);
434 }
435
436 CoordRepType a = p[1].SquaredEuclideanDistanceTo(p[2]);
437 CoordRepType b = p[0].SquaredEuclideanDistanceTo(p[2]);
438 CoordRepType c = p[1].SquaredEuclideanDistanceTo(p[0]);
439
440 weights[0] = a * ( b + c - a );
441 weights[1] = b * ( c + a - b );
442 weights[2] = c * ( a + b - c );
443
444 CoordRepType sum_weights = weights[0] + weights[1] + weights[2];
445
446 if ( sum_weights != 0. )
447 {
448 PointType oP;
449 oP.Fill(0.);
450
451 for ( i = 0; i < 3; i++ )
452 {
453 oP += p[i].GetVectorFromOrigin() * weights[i] / sum_weights;
454 }
455
456 return oP;
457 }
458 else
459 {
460 return p[0];
461 }
462 }
463
464 /** Evaluate the position of a given point inside the cell */
465 template< typename TCellInterface >
466 bool
467 TriangleCell< TCellInterface >
EvaluatePosition(CoordRepType * x,PointsContainer * points,CoordRepType * closestPoint,CoordRepType pcoord[3],double * minDist2,InterpolationWeightType * weights)468 ::EvaluatePosition(CoordRepType *x,
469 PointsContainer *points,
470 CoordRepType *closestPoint,
471 CoordRepType pcoord[3],
472 double *minDist2,
473 InterpolationWeightType *weights)
474 {
475 unsigned int i;
476 double dist2Point;
477 double dist2Line1;
478 double dist2Line2;
479 PointType closest;
480 PointType closestPoint1;
481 PointType closestPoint2;
482 PointType X(x);
483
484 if ( !points )
485 {
486 return false;
487 }
488
489 //
490 // Get the vertexes of this triangle
491 //
492 PointType pt1 = points->GetElement(m_PointIds[0]);
493 PointType pt2 = points->GetElement(m_PointIds[1]);
494 PointType pt3 = points->GetElement(m_PointIds[2]);
495
496 //
497 // Compute Vectors along the edges.
498 // These two vectors form a vector base for the 2D space of the triangle cell.
499 //
500 VectorType v12 = pt1 - pt2;
501 VectorType v32 = pt3 - pt2;
502
503 //
504 // Compute Vectors in the dual vector base inside the 2D space of the triangle
505 // cell.
506 // u12 is orthogonal to v32
507 // u32 is orthogonal to v12
508 //
509 const double dotproduct = v12 * v32;
510 VectorType u12 = v12 - v32 * ( dotproduct / v32.GetSquaredNorm() );
511 VectorType u32 = v32 - v12 * ( dotproduct / v12.GetSquaredNorm() );
512
513 //
514 // Add normalizations for making {u12,u32} a vector basis orthonormal to {v12,
515 // v32}.
516 //
517 u12 /= ( u12 * v12 );
518 u32 /= ( u32 * v32 );
519
520 //
521 // Project point to plane, by using the dual vector base
522 //
523 // Compute components of the input point in the 2D
524 // space defined by v12 and v32
525 //
526 VectorType xo = X - pt2;
527
528 const double u12p = xo * u12;
529 const double u32p = xo * u32;
530
531 VectorType x12 = v12 * u12p;
532 VectorType x32 = v32 * u32p;
533
534 //
535 // The projection of point X in the plane is cp
536 //
537 PointType cp = pt2 + x12 + x32;
538
539 //
540 // Compute barycentric coordinates in the Triangle
541 //
542 const double b1 = u12p;
543 const double b2 = 1.0 - u12p - u32p;
544 const double b3 = u32p;
545
546 //
547 // Test if the projected point is inside the cell.
548 //
549 // Zero with epsilon
550 const double zwe = -NumericTraits< double >::min();
551
552 //
553 // Since the three barycentric coordinates are interdependent
554 // only three tests should be necessary. That is, we only need
555 // to test against the equations of three lines (half-spaces).
556 //
557 if ( ( b1 >= zwe ) && ( b2 >= zwe ) && ( b3 >= zwe ) )
558 {
559 //
560 // This is the case when the point is inside the triangle
561 //projection distance
562 if ( closestPoint )
563 { // Compute the Distance 2 Between Points
564 *minDist2 = 0;
565 for ( i = 0; i < PointDimension; i++ )
566 {
567 const double val = cp[i] - x[i];
568 *minDist2 += val * val;
569 closestPoint[i] = cp[i];
570 }
571 }
572
573 if ( pcoord )
574 {
575 pcoord[0] = b1;
576 pcoord[1] = b2;
577 pcoord[2] = b3;
578 }
579
580 if ( weights )
581 {
582 weights[0] = b1;
583 weights[1] = b2;
584 weights[2] = b3;
585 }
586
587 return true;
588 }
589 else
590 {
591 if ( closestPoint )
592 {
593 double lt; // parameter along the line (not used)
594 if ( b1 < 0.0 && b2 < 0.0 )
595 {
596 dist2Point = 0;
597 for ( i = 0; i < PointDimension; i++ )
598 {
599 dist2Point += (x[i] - pt3[i]) * (x[i] - pt3[i]);
600 }
601 dist2Line1 = this->DistanceToLine(x, pt1, pt3, lt, closestPoint1);
602 dist2Line2 = this->DistanceToLine(x, pt3, pt2, lt, closestPoint2);
603 if ( dist2Point < dist2Line1 )
604 {
605 *minDist2 = dist2Point;
606 closest = pt3;
607 }
608 else
609 {
610 *minDist2 = dist2Line1;
611 closest = closestPoint1;
612 }
613 if ( dist2Line2 < *minDist2 )
614 {
615 *minDist2 = dist2Line2;
616 closest = closestPoint2;
617 }
618 for ( i = 0; i < PointDimension; i++ )
619 {
620 closestPoint[i] = closest[i];
621 }
622 for(; i < 3; i++ )
623 {
624 closestPoint[i] = 0.;
625 }
626 }
627 else if ( b2 < 0.0 && b3 < 0.0 )
628 {
629 dist2Point = 0;
630 for ( i = 0; i < PointDimension; i++ )
631 {
632 dist2Point += (x[i] - pt1[i]) * (x[i] - pt1[i]);
633 }
634 dist2Line1 = this->DistanceToLine(x, pt1, pt3, lt, closestPoint1);
635 dist2Line2 = this->DistanceToLine(x, pt1, pt2, lt, closestPoint2);
636 if ( dist2Point < dist2Line1 )
637 {
638 *minDist2 = dist2Point;
639 closest = pt1;
640 }
641 else
642 {
643 *minDist2 = dist2Line1;
644 closest = closestPoint1;
645 }
646 if ( dist2Line2 < *minDist2 )
647 {
648 *minDist2 = dist2Line2;
649 closest = closestPoint2;
650 }
651 for ( i = 0; i < PointDimension; i++ )
652 {
653 closestPoint[i] = closest[i];
654 }
655 for(; i < 3; i++ )
656 {
657 closestPoint[i] = 0.;
658 }
659 }
660 else if ( b1 < 0.0 && b3 < 0.0 )
661 {
662 dist2Point = 0;
663 for ( i = 0; i < PointDimension; i++ )
664 {
665 dist2Point += ( x[i] - pt2[i] ) * ( x[i] - pt2[i] );
666 }
667 dist2Line1 = this->DistanceToLine(x, pt2, pt3, lt, closestPoint1);
668 dist2Line2 = this->DistanceToLine(x, pt1, pt2, lt, closestPoint2);
669 if ( dist2Point < dist2Line1 )
670 {
671 *minDist2 = dist2Point;
672 closest = pt2;
673 }
674 else
675 {
676 *minDist2 = dist2Line1;
677 closest = closestPoint1;
678 }
679 if ( dist2Line2 < *minDist2 )
680 {
681 *minDist2 = dist2Line2;
682 closest = closestPoint2;
683 }
684 for ( i = 0; i < PointDimension; i++ )
685 {
686 closestPoint[i] = closest[i];
687 }
688 for(; i < 3; i++ )
689 {
690 closestPoint[i] = 0.;
691 }
692 }
693 else if ( b1 < 0.0 )
694 {
695 *minDist2 = this->DistanceToLine(x, pt2, pt3, lt, closestPoint);
696 }
697 else if ( b2 < 0.0 )
698 {
699 *minDist2 = this->DistanceToLine(x, pt1, pt3, lt, closestPoint);
700 }
701 else if ( b3 < 0.0 )
702 {
703 *minDist2 = this->DistanceToLine(x, pt1, pt2, lt, closestPoint);
704 }
705 }
706 if ( pcoord )
707 {
708 pcoord[0] = b1;
709 pcoord[1] = b2;
710 pcoord[2] = b3;
711 }
712 //Just fall through to default return false;
713 }
714 return false; //Default case that should never be reached.
715 }
716 } // end namespace itk
717
718 #endif
719