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