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