KWStyle - itkTriangleCell.txx

Matrix View

Description

---

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

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:51PM

© Kitware Inc.