1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkCubicLine.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkCubicLine.h"
16
17 #include "vtkNonLinearCell.h"
18 #include "vtkCellArray.h"
19 #include "vtkCellData.h"
20 #include "vtkCell.h"
21 #include "vtkMath.h"
22 #include "vtkObjectFactory.h"
23 #include "vtkIncrementalPointLocator.h"
24 #include "vtkLine.h"
25 #include "vtkDoubleArray.h"
26 #include "vtkPoints.h"
27 #include "vtkPointData.h"
28
29 vtkStandardNewMacro(vtkCubicLine);
30
31 //----------------------------------------------------------------------------
32 // Construct the line with four points.
vtkCubicLine()33 vtkCubicLine::vtkCubicLine()
34 {
35 this->Scalars = vtkDoubleArray::New();
36 this->Scalars->SetNumberOfTuples(4);
37 this->Points->SetNumberOfPoints(4);
38 this->PointIds->SetNumberOfIds(4);
39 for (int i = 0; i < 4; i++)
40 {
41 this->Points->SetPoint(i, 0.0, 0.0, 0.0);
42 this->PointIds->SetId(i,0);
43 }
44 this->Line = vtkLine::New();
45 }
46 //----------------------------------------------------------------------------
47 // Delete the Line
~vtkCubicLine()48 vtkCubicLine::~vtkCubicLine()
49 {
50 this->Line->Delete();
51 this->Scalars->Delete();
52 }
53
54
55 //----------------------------------------------------------------------------
EvaluatePosition(const double x[3],double closestPoint[3],int & subId,double pcoords[3],double & minDist2,double weights[])56 int vtkCubicLine::EvaluatePosition(const double x[3], double closestPoint[3],
57 int& subId, double pcoords[3],
58 double& minDist2, double weights[])
59 {
60 double closest[3];
61 double pc[3], dist2;
62 int ignoreId, i, returnStatus, status;
63 double lineWeights[2];
64
65 pcoords[1] = pcoords[2] = 0.0;
66
67 returnStatus = -1;
68 weights[0] = 0.0;
69 for (minDist2=VTK_DOUBLE_MAX,i=0; i < 3; i++)
70 {
71 if ( i == 0)
72 {
73 this->Line->Points->SetPoint(0,this->Points->GetPoint(0));
74 this->Line->Points->SetPoint(1,this->Points->GetPoint(2));
75 }
76 else if (i == 1)
77 {
78 this->Line->Points->SetPoint(0,this->Points->GetPoint(2));
79 this->Line->Points->SetPoint(1,this->Points->GetPoint(3));
80 }
81 else
82 {
83 this->Line->Points->SetPoint(0,this->Points->GetPoint(3));
84 this->Line->Points->SetPoint(1,this->Points->GetPoint(1));
85 }
86
87
88 status = this->Line->EvaluatePosition(x,closest,ignoreId,pc,
89 dist2,lineWeights);
90 if ( status != -1 && dist2 < minDist2 )
91 {
92 returnStatus = status;
93 minDist2 = dist2;
94 subId = i;
95 pcoords[0] = pc[0];
96 }
97 }
98
99 // adjust parametric coordinate
100 if ( returnStatus != -1 )
101 {
102 if ( subId == 0 ) //first part : -1 <= pcoords <= -1/3
103 {
104 pcoords[0] = pcoords[0]*(2.0/3.0) - 1;
105 }
106 else if ( subId == 1 ) //second part : -1/3 <= pcoords <= 1/3
107 {
108 pcoords[0] = pcoords[0]*(2.0/3.0) -(1.0/3.0) ;
109 }
110 else // third part : 1/3 <= pcoords <= 1
111 {
112 pcoords[0] = pcoords[0]*(2.0/3.0) + (1.0/3.0);
113 }
114 if(closestPoint!=nullptr)
115 {
116 // Compute both closestPoint and weights
117 this->EvaluateLocation(subId,pcoords,closestPoint,weights);
118 }
119 else
120 {
121 // Compute weights only
122 this->InterpolationFunctions(pcoords,weights);
123 }
124 }
125
126 return returnStatus;
127 }
128
129
130
131 //----------------------------------------------------------------------------
EvaluateLocation(int & vtkNotUsed (subId),const double pcoords[3],double x[3],double * weights)132 void vtkCubicLine::EvaluateLocation(int& vtkNotUsed(subId), const double pcoords[3],
133 double x[3], double *weights)
134 {
135 int i;
136 double a0[3], a1[3], a2[3], a3[3];
137 this->Points->GetPoint(0, a0);
138 this->Points->GetPoint(1, a1);
139 this->Points->GetPoint(2, a2); //first midside node
140 this->Points->GetPoint(3, a3); //second midside node
141
142 this->InterpolationFunctions(pcoords,weights);
143
144
145 for (i=0; i<3; i++)
146 {
147 x[i] = a0[i]*weights[0] + a1[i]*weights[1] + a2[i]*weights[2] + a3[i]*weights[3];
148 }
149
150 }
151
152
153
154
155
156 //----------------------------------------------------------------------------
CellBoundary(int vtkNotUsed (subId),const double pcoords[3],vtkIdList * pts)157 int vtkCubicLine::CellBoundary(int vtkNotUsed(subId), const double pcoords[3],
158 vtkIdList *pts)
159 {
160 pts->SetNumberOfIds(1);
161
162 if ( pcoords[0] >= 0.0 )
163 {
164 pts->SetId(0,this->PointIds->GetId(1)); // The edge points IDs are 0 and 1.
165 if ( pcoords[0] > 1.0 )
166 {
167 return 0;
168 }
169 else
170 {
171 return 1;
172 }
173 }
174 else
175 {
176 pts->SetId(0,this->PointIds->GetId(0));
177 if ( pcoords[0] < -1.0 )
178 {
179 return 0;
180 }
181 else
182 {
183 return 1;
184 }
185 }
186 }
187
188
189 //LinearLines for the Contour and the Clip Algorithm
190 //-----------------------------------------------------------------------------
191 static int LinearLines[3][2] = { {0,2}, {2,3}, {3,1} };
192
193
Contour(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * verts,vtkCellArray * lines,vtkCellArray * polys,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd)194 void vtkCubicLine::Contour(double value, vtkDataArray *cellScalars,
195 vtkIncrementalPointLocator *locator, vtkCellArray *verts,
196 vtkCellArray *lines, vtkCellArray *polys,
197 vtkPointData *inPd, vtkPointData *outPd,
198 vtkCellData *inCd, vtkIdType cellId,
199 vtkCellData *outCd)
200 {
201 for (int i=0; i < 3; i++) //for each subdivided line
202 {
203 for ( int j=0; j<2; j++) //for each of the four vertices of the line
204 {
205 this->Line->Points->SetPoint(j,this->Points->GetPoint(LinearLines[i][j]));
206 this->Line->PointIds->SetId(j,this->PointIds->GetId(LinearLines[i][j]));
207 this->Scalars->SetValue(j,cellScalars->GetTuple1(LinearLines[i][j]));
208 }
209 this->Line->Contour(value, this->Scalars, locator, verts,
210 lines, polys, inPd, outPd, inCd, cellId, outCd);
211 }
212 }
213
214
215
216
217 //----------------------------------------------------------------------------
218 // Line-line intersection. Intersection has to occur within [0,1] parametric
219 // coordinates and with specified tolerance.
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId)220 int vtkCubicLine::IntersectWithLine(const double p1[3], const double p2[3], double tol, double& t,
221 double x[3], double pcoords[3], int& subId)
222 {
223
224
225 int subTest, numLines=3;
226
227 for (subId=0; subId < numLines; subId++)
228 {
229 if ( subId == 0)
230 {
231 this->Line->Points->SetPoint(0,this->Points->GetPoint(0));
232 this->Line->Points->SetPoint(1,this->Points->GetPoint(2));
233 }
234 else if (subId == 1)
235 {
236 this->Line->Points->SetPoint(0,this->Points->GetPoint(2));
237 this->Line->Points->SetPoint(1,this->Points->GetPoint(3));
238 }
239 else
240 {
241 this->Line->Points->SetPoint(0,this->Points->GetPoint(3));
242 this->Line->Points->SetPoint(1,this->Points->GetPoint(1));
243 }
244
245 if ( this->Line->IntersectWithLine(p1, p2, tol, t, x, pcoords, subTest) )
246 {
247 // adjust parametric coordinate
248
249 if ( subId == 0 ) //first part : -1 <= pcoords <= -1/3
250 {
251 pcoords[0] = pcoords[0]*(2.0/3.0) - 1;
252 }
253 else if ( subId == 1 ) //second part : -1/3 <= pcoords <= 1/3
254 {
255 pcoords[0] = pcoords[0]*(2.0/3.0) - (1.0/3.0) ;
256 }
257 else // third part : 1/3 <= pcoords <= 1
258 {
259 pcoords[0] = pcoords[0]*(2.0/3.0) + (1.0/3.0);
260 }
261
262 return 1;
263 }
264 }
265
266 return 0;
267 }
268
269
270
271
272 //----------------------------------------------------------------------------
Triangulate(int vtkNotUsed (index),vtkIdList * ptIds,vtkPoints * pts)273 int vtkCubicLine::Triangulate(int vtkNotUsed(index), vtkIdList *ptIds,
274 vtkPoints *pts)
275 {
276 pts->Reset();
277 ptIds->Reset();
278
279 // The first line
280 ptIds->InsertId(0,this->PointIds->GetId(0));
281 pts->InsertPoint(0,this->Points->GetPoint(0));
282
283 ptIds->InsertId(1,this->PointIds->GetId(2));
284 pts->InsertPoint(1,this->Points->GetPoint(2));
285
286 // The second line
287 ptIds->InsertId(2,this->PointIds->GetId(2));
288 pts->InsertPoint(2,this->Points->GetPoint(2));
289
290 ptIds->InsertId(3,this->PointIds->GetId(3));
291 pts->InsertPoint(3,this->Points->GetPoint(3));
292
293
294 // The third line
295 ptIds->InsertId(4,this->PointIds->GetId(3));
296 pts->InsertPoint(4,this->Points->GetPoint(3));
297
298 ptIds->InsertId(5,this->PointIds->GetId(1));
299 pts->InsertPoint(5,this->Points->GetPoint(1));
300 return 1;
301 }
302
303
304 //----------------------------------------------------------------------------
Derivatives(int vtkNotUsed (subId),const double pcoords[3],const double * values,int dim,double * derivs)305 void vtkCubicLine::Derivatives(int vtkNotUsed(subId),
306 const double pcoords[3],
307 const double *values,
308 int dim,
309 double *derivs)
310 {
311 double v0, v1, v2, v3; // Local coordinates of Each point.
312 double v10[3], lenX; // Reesentation of local Axis
313 double x0[3], x1[3], x2[3], x3[3]; // Points of the model
314 double vec20[3], vec30[3]; // Normal and vector of each point.
315 double J; // Jacobian Matrix
316 double JI; // Inverse of the Jacobian Matrix
317 double funcDerivs[4], sum, dBydx; // Derivated values
318
319
320 // Project points of vtkCubicLine into a 1D system
321 this->Points->GetPoint(0, x0);
322 this->Points->GetPoint(1, x1);
323 this->Points->GetPoint(2, x2);
324 this->Points->GetPoint(3, x3);
325
326 for (int i=0; i < 3; i++) // Compute the vector for each point
327 {
328 v10[i] = x1[i] - x0[i];
329 vec20[i] = x2[i] - x0[i];
330 vec30[i] = x3[i] - x0[i];
331 }
332
333
334 if ( (lenX=vtkMath::Normalize(v10)) <= 0.0 ) //degenerate
335 {
336 for (int j=0; j < dim; j++ )
337 {
338 for (int i=0; i < 3; i++ )
339 {
340 derivs[j*dim + i] = 0.0;
341 }
342 }
343 return;
344 }
345
346 v0 = 0.0; //convert points to 1D (i.e., local system)
347
348 v1 = lenX;
349
350 v2 = vtkMath::Dot(vec20,v10);
351
352 v3 = vtkMath::Dot(vec30,v10);
353
354 this->InterpolationDerivs(pcoords, funcDerivs);
355
356 J = v0*funcDerivs[0] + v1*funcDerivs[1] +
357 v2*funcDerivs[2] + v3*funcDerivs[3];
358
359 // Compute inverse Jacobian, return if Jacobian is singular
360
361 if(J != 0)
362 {
363 JI = 1.0 / J;
364 }
365 else
366 {
367 for (int j=0; j < dim; j++ )
368 {
369 for (int i=0; i < 3; i++ )
370 {
371 derivs[j*dim + i] = 0.0;
372 }
373 }
374 return;
375 }
376
377 // Loop over "dim" derivative values. For each set of values,
378 // compute derivatives
379 // in local system and then transform into modelling system.
380 // First compute derivatives in local x' coordinate system
381 for (int j=0; j < dim; j++ )
382 {
383 sum = 0.0;
384 for (int i=0; i < 4; i++) //loop over interp. function derivatives
385 {
386 sum += funcDerivs[i] * values[dim*i + j];
387 }
388
389 dBydx = sum*JI;
390
391 // Transform into global system (dot product with global axes)
392 derivs[3*j] = dBydx * v10[0] ;
393 derivs[3*j + 1] = dBydx * v10[1] ;
394 derivs[3*j + 2] = dBydx * v10[2] ;
395 }
396 }
397
398
399
400
401 //--------------------------------------------------------------------------
402 // Clip this line using scalar value provided. Like contouring, except
403 // that it cuts the line to produce other lines.
Clip(double value,vtkDataArray * cellScalars,vtkIncrementalPointLocator * locator,vtkCellArray * lines,vtkPointData * inPd,vtkPointData * outPd,vtkCellData * inCd,vtkIdType cellId,vtkCellData * outCd,int insideOut)404 void vtkCubicLine::Clip(double value, vtkDataArray *cellScalars,
405 vtkIncrementalPointLocator *locator, vtkCellArray *lines,
406 vtkPointData *inPd, vtkPointData *outPd,
407 vtkCellData *inCd, vtkIdType cellId, vtkCellData *outCd,
408 int insideOut)
409 {
410 for (int i=0; i < 3; i++) //for each subdivided line
411 {
412 for ( int j=0; j<2; j++) //for each of the four vertices of the line
413 {
414 this->Line->Points->SetPoint(j,this->Points->GetPoint(LinearLines[i][j]));
415 this->Line->PointIds->SetId(j,this->PointIds->GetId(LinearLines[i][j]));
416 this->Scalars->SetValue(j,cellScalars->GetTuple1(LinearLines[i][j]));
417 }
418 this->Line->Clip(value, this->Scalars, locator, lines, inPd, outPd,
419 inCd, cellId, outCd, insideOut);
420 }
421 }
422
423
424 //----------------------------------------------------------------------------
425 //
426 // Compute interpolation functions
427 //
InterpolationFunctions(const double pcoords[3],double weights[4])428 void vtkCubicLine::InterpolationFunctions(const double pcoords[3], double weights[4]) // N2 and N3 are the middle points
429 {
430 // pcoords[0] = t, weights need to be set in accordance with the definition of the standard cubic line finite element
431 double t = pcoords[0];
432
433 weights[0] = (9.0/16.0)*(1.0 - t)*(t + (1.0/3.0))*(t - (1.0/3.0));
434 weights[1] = (-9.0/16.0)*(1.0 + t)*((1.0/3.0) - t)*(t + (1.0/3.0));
435 weights[2] = (27.0/16.0)*(t - 1.0)*(t + 1.0)*(t - (1.0/3.0));
436 weights[3] = (-27.0/16.0)*(t - 1.0)*(t + 1.0)*(t + (1.0/3.0));
437 }
438
439
440
441
442 //----------------------------------------------------------------------------
InterpolationDerivs(const double pcoords[3],double derivs[4])443 void vtkCubicLine::InterpolationDerivs(const double pcoords[3], double derivs[4]) //N2 and N3 are the middle points
444 {
445 double t = pcoords[0];
446
447 derivs[0] = (1.0/16.0)*(1.0 + 18.0*t - 27.0*t*t);
448 derivs[1] = (1.0/16.0)*(-1.0 + 18.0*t + 27.0*t*t);
449 derivs[2] = (1.0/16.0)*(-27.0 - 18.0*t + 81.0*t*t);
450 derivs[3] = (1.0/16.0)*(27.0 - 18.0*t - 81.0*t*t);
451
452 }
453
454
455
456 //----------------------------------------------------------------------------
457 static double vtkCubicLineCellPCoords[12] = {-1.0,0.0,0.0, 1.0,0.0,0.0, -(1.0/3.0),0.0,0.0, (1.0/3.0),0.0,0.0};
GetParametricCoords()458 double *vtkCubicLine::GetParametricCoords()
459 {
460 return vtkCubicLineCellPCoords;
461 }
462
463
464 //----------------------------------------------------------------------------
GetParametricDistance(const double pcoords[3])465 double vtkCubicLine::GetParametricDistance(const double pcoords[3])
466 {
467
468 double pc;
469
470 pc = pcoords[0];
471
472 if ( pc <= -1.0)
473 {
474 return pc * (-1.0) - 1.0;
475 }
476 else if (pc >= 1.0)
477 {
478 return pc - 1.0;
479 }
480
481 return pc; // the parametric coordintate lies between -1.0 and 1.0.
482 }
483
484
485 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)486 void vtkCubicLine::PrintSelf(ostream& os, vtkIndent indent)
487 {
488 this->Superclass::PrintSelf(os,indent);
489 os << indent << "Line: " << this->Line << endl;
490 }
491