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 itkQuadrilateralCell_hxx
19 #define itkQuadrilateralCell_hxx
20 #include "itkQuadrilateralCell.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 QuadrilateralCell< TCellInterface >
MakeCopy(CellAutoPointer & cellPointer) const32 ::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 QuadrilateralCell< TCellInterface >
GetDimension() const45 ::GetDimension() 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 QuadrilateralCell< TCellInterface >
GetNumberOfPoints() const57 ::GetNumberOfPoints() 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 QuadrilateralCell< TCellInterface >::CellFeatureCount
68 QuadrilateralCell< TCellInterface >
GetNumberOfBoundaryFeatures(int dimension) const69 ::GetNumberOfBoundaryFeatures(int dimension) const
70 {
71 switch ( dimension )
72 {
73 case 0:
74 return GetNumberOfVertices();
75 case 1:
76 return GetNumberOfEdges();
77 default:
78 return 0;
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 QuadrilateralCell< TCellInterface >
GetBoundaryFeature(int dimension,CellFeatureIdentifier featureId,CellAutoPointer & cellPointer)91 ::GetBoundaryFeature(int dimension, CellFeatureIdentifier featureId, 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 QuadrilateralCell< TCellInterface >
SetPointIds(PointIdConstIterator first)131 ::SetPointIds(PointIdConstIterator first)
132 {
133 PointIdConstIterator ii(first);
134
135 for ( unsigned int i = 0; i < Self::NumberOfPoints; ++i )
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 QuadrilateralCell< TCellInterface >
SetPointIds(PointIdConstIterator first,PointIdConstIterator last)151 ::SetPointIds(PointIdConstIterator first, PointIdConstIterator last)
152 {
153 int localId = 0;
154 PointIdConstIterator ii(first);
155
156 while ( ii != last )
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 QuadrilateralCell< 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 QuadrilateralCell< TCellInterface >::PointIdIterator
180 QuadrilateralCell< 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 QuadrilateralCell< TCellInterface >::PointIdConstIterator
193 QuadrilateralCell< 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 QuadrilateralCell< TCellInterface >::PointIdIterator
205 QuadrilateralCell< 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 QuadrilateralCell< TCellInterface >::PointIdConstIterator
218 QuadrilateralCell< TCellInterface >
PointIdsEnd() const219 ::PointIdsEnd() const
220 {
221 return &m_PointIds[Self::NumberOfPoints - 1] + 1;
222 }
223
224 /**
225 * Quadrilateral-specific:
226 * Get the number of vertices defining the quadrilateral.
227 */
228 template< typename TCellInterface >
229 typename QuadrilateralCell< TCellInterface >::CellFeatureCount
230 QuadrilateralCell< TCellInterface >
GetNumberOfVertices() const231 ::GetNumberOfVertices() const
232 {
233 return NumberOfVertices;
234 }
235
236 /**
237 * Quadrilateral-specific:
238 * Get the number of edges defined for the quadrilateral.
239 */
240 template< typename TCellInterface >
241 typename QuadrilateralCell< TCellInterface >::CellFeatureCount
242 QuadrilateralCell< TCellInterface >
GetNumberOfEdges() const243 ::GetNumberOfEdges() const
244 {
245 return Self::NumberOfEdges;
246 }
247
248 /**
249 * Quadrilateral-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 QuadrilateralCell< 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 * Quadrilateral-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 QuadrilateralCell< 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 /** Evaluate the position inside the cell */
286 template< typename TCellInterface >
287 bool
288 QuadrilateralCell< TCellInterface >
EvaluatePosition(CoordRepType * x,PointsContainer * points,CoordRepType * closestPoint,CoordRepType pcoord[CellDimension],double * dist2,InterpolationWeightType * weight)289 ::EvaluatePosition(CoordRepType *x,
290 PointsContainer *points,
291 CoordRepType *closestPoint,
292 CoordRepType pcoord[CellDimension],
293 double *dist2,
294 InterpolationWeightType *weight)
295 {
296 static constexpr int ITK_QUAD_MAX_ITERATION = 10;
297 static constexpr double ITK_QUAD_CONVERGED = 1.e-03;
298 static constexpr double ITK_DIVERGED = 1.e6;
299
300 int iteration, converged;
301 double params[CellDimension];
302 double fcol[CellDimension];
303 double rcol[CellDimension];
304 double scol[CellDimension];
305 double d;
306 PointType pt;
307 CoordRepType derivs[NumberOfDerivatives];
308 InterpolationWeightType weights[NumberOfPoints];
309
310 // set initial position for Newton's method
311 int subId = 0;
312 CoordRepType pcoords[CellDimension];
313
314 pcoords[0] = pcoords[1] = params[0] = params[1] = 0.5;
315
316 // NOTE: Point x is here assumed to lie on the plane of Quad. Otherwise, (FIXME)
317 // - Get normal for quadrilateral, using its 3 corners
318 // - Project point x onto Quad plane using this normal
319 // See vtkQuad for this: ComputeNormal (this, pt1, pt2, pt3, n); vtkPlane::ProjectPoint(x,pt1,n,cp);
320
321 // enter iteration loop
322 for ( iteration = converged = 0;
323 !converged && ( iteration < ITK_QUAD_MAX_ITERATION ); iteration++ )
324 {
325 // calculate element interpolation functions and derivatives
326 this->InterpolationFunctions(pcoords, weights);
327 this->InterpolationDerivs(pcoords, derivs);
328
329 // calculate newton functions
330 for ( unsigned int i = 0; i < CellDimension; ++i )
331 {
332 fcol[i] = rcol[i] = scol[i] = 0.0;
333 }
334 for ( unsigned int i = 0; i < NumberOfPoints; ++i )
335 {
336 pt = points->GetElement(m_PointIds[i]);
337 // using the projection normal n, one can choose which 2 axes to use out of 3
338 // any 2 should work, so (not having n) we use [x,y] (also assuming 2D use of QuadCell)
339 // if we compute n, one can use the closest two indices as in vtkQuad
340 for ( unsigned int j = 0; j < CellDimension; ++j )
341 {
342 fcol[j] += pt[j] * weights[i];
343 rcol[j] += pt[j] * derivs[i];
344 scol[j] += pt[j] * derivs[i + NumberOfPoints];
345 }
346 }
347
348 for ( unsigned int i = 0; i < CellDimension; ++i )
349 {
350 fcol[i] -= x[i];
351 }
352
353 // compute determinants and generate improvements
354 vnl_matrix_fixed< CoordRepType, CellDimension, CellDimension > mat;
355 for ( unsigned int i = 0; i < CellDimension; ++i )
356 {
357 mat.put(0, i, rcol[i]);
358 mat.put(1, i, scol[i]);
359 }
360
361 d = vnl_determinant(mat);
362 //d=vtkMath::Determinant2x2(rcol,scol);
363 if ( std::abs(d) < 1.e-20 )
364 {
365 return false;
366 }
367
368 vnl_matrix_fixed< CoordRepType, CellDimension, CellDimension > mat1;
369 for ( unsigned int i = 0; i < CellDimension; ++i )
370 {
371 mat1.put(0, i, fcol[i]);
372 mat1.put(1, i, scol[i]);
373 }
374
375 vnl_matrix_fixed< CoordRepType, CellDimension, CellDimension > mat2;
376 for ( unsigned int i = 0; i < CellDimension; ++i )
377 {
378 mat2.put(0, i, rcol[i]);
379 mat2.put(1, i, fcol[i]);
380 }
381
382 pcoords[0] = params[0] - vnl_determinant(mat1) / d;
383 pcoords[1] = params[1] - vnl_determinant(mat2) / d;
384
385 if ( pcoord )
386 {
387 pcoord[0] = pcoords[0];
388 pcoord[1] = pcoords[1];
389 }
390
391 // check for convergence
392 if ( ( ( std::abs(pcoords[0] - params[0]) ) < ITK_QUAD_CONVERGED )
393 && ( ( std::abs(pcoords[1] - params[1]) ) < ITK_QUAD_CONVERGED ) )
394 {
395 converged = 1;
396 }
397
398 // Test for bad divergence (S.Hirschberg 11.12.2001)
399 else if ( ( std::abs(pcoords[0]) > ITK_DIVERGED )
400 || ( std::abs(pcoords[1]) > ITK_DIVERGED ) )
401 {
402 return -1;
403 }
404
405 // if not converged, repeat
406 else
407 {
408 params[0] = pcoords[0];
409 params[1] = pcoords[1];
410 }
411 }
412
413 // if not converged, set the parametric coordinates to arbitrary values
414 // outside of element
415 if ( !converged )
416 {
417 return false;
418 }
419
420 this->InterpolationFunctions(pcoords, weights);
421
422 if ( weight )
423 {
424 for ( unsigned int i = 0; i < NumberOfPoints; ++i )
425 {
426 weight[i] = weights[i];
427 }
428 }
429
430 if ( pcoords[0] >= -0.001 && pcoords[0] <= 1.001
431 && pcoords[1] >= -0.001 && pcoords[1] <= 1.001 )
432 {
433 if ( closestPoint )
434 {
435 closestPoint[0] = x[0];
436 closestPoint[1] = x[1];
437 *dist2 = 0.0; //inside quadrilateral
438 }
439 return true;
440 }
441 else
442 {
443 CoordRepType pc[CellDimension], w[NumberOfPoints];
444 if ( closestPoint )
445 {
446 for ( unsigned int i = 0; i < CellDimension; ++i ) //only approximate ??
447 {
448 if ( pcoords[i] < 0.0 )
449 {
450 pc[i] = 0.0;
451 }
452 else if ( pcoords[i] > 1.0 )
453 {
454 pc[i] = 1.0;
455 }
456 else
457 {
458 pc[i] = pcoords[i];
459 }
460 }
461 this->EvaluateLocation(subId, points, pc, closestPoint, (InterpolationWeightType *)w);
462
463 *dist2 = 0;
464 for ( unsigned int i = 0; i < CellDimension; ++i )
465 {
466 *dist2 += ( closestPoint[i] - x[i] ) * ( closestPoint[i] - x[i] );
467 }
468 }
469 return false;
470 }
471 }
472
473 /** Compute iso-parametric interpolation functions */
474 template< typename TCellInterface >
475 void
476 QuadrilateralCell< TCellInterface >
InterpolationFunctions(const CoordRepType pointCoords[CellDimension],InterpolationWeightType weights[NumberOfPoints])477 ::InterpolationFunctions(const CoordRepType pointCoords[CellDimension], InterpolationWeightType weights[NumberOfPoints])
478 {
479 const double rm = 1. - pointCoords[0];
480 const double sm = 1. - pointCoords[1];
481
482 weights[0] = rm * sm;
483 weights[1] = pointCoords[0] * sm;
484 weights[2] = pointCoords[0] * pointCoords[1];
485 weights[3] = rm * pointCoords[1];
486 }
487
488 /** Compute iso-parametric interpolation functions */
489 template< typename TCellInterface >
490 void
491 QuadrilateralCell< TCellInterface >
InterpolationDerivs(const CoordRepType pointCoords[CellDimension],CoordRepType derivs[NumberOfDerivatives])492 ::InterpolationDerivs(const CoordRepType pointCoords[CellDimension], CoordRepType derivs[NumberOfDerivatives])
493 {
494 const double rm = 1. - pointCoords[0];
495 const double sm = 1. - pointCoords[1];
496
497 // r-derivatives
498 derivs[0] = -sm;
499 derivs[1] = sm;
500 derivs[2] = pointCoords[1];
501 derivs[3] = -pointCoords[1];
502 // s-derivatives
503 derivs[4] = -rm;
504 derivs[5] = -pointCoords[0];
505 derivs[6] = pointCoords[0];
506 derivs[7] = rm;
507 }
508
509 /** Evaluate the location inside the cell */
510 template< typename TCellInterface >
511 void
512 QuadrilateralCell< TCellInterface >
EvaluateLocation(int & itkNotUsed (subId),const PointsContainer * points,const CoordRepType pointCoords[PointDimension],CoordRepType x[PointDimension],InterpolationWeightType * weights)513 ::EvaluateLocation(int & itkNotUsed(subId), const PointsContainer *points, const CoordRepType pointCoords[PointDimension],
514 CoordRepType x[PointDimension], InterpolationWeightType *weights)
515 {
516 this->InterpolationFunctions(pointCoords, weights);
517
518 for ( unsigned int ii = 0; ii < PointDimension; ++ii )
519 {
520 x[ii] = NumericTraits< CoordRepType >::ZeroValue();
521 }
522
523 for ( unsigned int ii = 0; ii < NumberOfPoints; ++ii )
524 {
525 const PointType & point = points->GetElement(m_PointIds[ii]);
526
527 for ( unsigned int jj = 0; jj < PointDimension; ++jj )
528 {
529 x[jj] += point[jj] * weights[ii];
530 }
531 }
532 }
533
534 } // end namespace itk
535
536 #endif
537