1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkSplineGraphEdges.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 /*-------------------------------------------------------------------------
16 Copyright 2008 Sandia Corporation.
17 Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
18 the U.S. Government retains certain rights in this software.
19 -------------------------------------------------------------------------*/
20 #include "vtkSplineGraphEdges.h"
21
22 #include "vtkCardinalSpline.h"
23 #include "vtkCommand.h"
24 #include "vtkGraph.h"
25 #include "vtkInformation.h"
26 #include "vtkInformationVector.h"
27 #include "vtkMath.h"
28 #include "vtkObjectFactory.h"
29
30 vtkStandardNewMacro(vtkSplineGraphEdges);
31 vtkCxxSetObjectMacro(vtkSplineGraphEdges, Spline, vtkSpline);
32
33 namespace {
34 // N-function defined at:
35 // http://mathworld.wolfram.com/B-Spline.html
36 // optimized for j = 3.
CubicSpline(vtkIdType i,double * k,double t)37 double CubicSpline(vtkIdType i, double* k, double t)
38 {
39 if (t >= k[i] && t < k[i+1])
40 {
41 double denom = (k[i+3]-k[i])*(k[i+2]-k[i])*(k[i+1]-k[i]);
42 if (denom == 0.0) return 0.0;
43 double temp = t - k[i];
44 return temp*temp*temp/denom;
45 }
46 if (t >= k[i+1] && t < k[i+2])
47 {
48 double denom1 = (k[i+3]-k[i])*(k[i+2]-k[i])*(k[i+2]-k[i+1]);
49 double term1;
50 if (denom1 == 0.0)
51 {
52 term1 = 0.0;
53 }
54 else
55 {
56 term1 = (t-k[i])*(t-k[i])*(k[i+2]-t)/denom1;
57 }
58
59 double denom2 = (k[i+3]-k[i])*(k[i+3]-k[i+1])*(k[i+2]-k[i+1]);
60 double term2;
61 if (denom2 == 0.0)
62 {
63 term2 = 0.0;
64 }
65 else
66 {
67 term2 = (t-k[i])*(k[i+3]-t)*(t-k[i+1])/denom2;
68 }
69
70 double denom3 = (k[i+4]-k[i+1])*(k[i+3]-k[i+1])*(k[i+2]-k[i+1]);
71 double term3;
72 if (denom3 == 0.0)
73 {
74 term3 = 0.0;
75 }
76 else
77 {
78 term3 = (k[i+4]-t)*(t-k[i+1])*(t-k[i+1])/denom3;
79 }
80
81 return term1 + term2 + term3;
82 }
83 if (t >= k[i+2] && t < k[i+3])
84 {
85 double denom1 = (k[i+3]-k[i])*(k[i+3]-k[i+1])*(k[i+3]-k[i+2]);
86 double term1;
87 if (denom1 == 0.0)
88 {
89 term1 = 0.0;
90 }
91 else
92 {
93 term1 = (t-k[i])*(k[i+3]-t)*(k[i+3]-t)/denom1;
94 }
95
96 double denom2 = (k[i+4]-k[i+1])*(k[i+3]-k[i+1])*(k[i+3]-k[i+2]);
97 double term2;
98 if (denom2 == 0.0)
99 {
100 term2 = 0.0;
101 }
102 else
103 {
104 term2 = (k[i+4]-t)*(t-k[i+1])*(k[i+3]-t)/denom2;
105 }
106
107 double denom3 = (k[i+4]-k[i+1])*(k[i+4]-k[i+2])*(k[i+3]-k[i+2]);
108 double term3;
109 if (denom3 == 0.0)
110 {
111 term3 = 0.0;
112 }
113 else
114 {
115 term3 = (k[i+4]-t)*(k[i+4]-t)*(t-k[i+2])/denom3;
116 }
117
118 return term1 + term2 + term3;
119 }
120 if (t >= k[i+3] && t < k[i+4])
121 {
122 double denom = (k[i+4]-k[i+1])*(k[i+4]-k[i+2])*(k[i+4]-k[i+3]);
123 if (denom == 0.0) return 0.0;
124 double temp = k[i+4] - t;
125 return temp*temp*temp/denom;
126 }
127 return 0.0;
128 }
129
130 // Slow, recursive version of N-function defined:
131 // http://mathworld.wolfram.com/B-Spline.html
132 #if 0
133 double N(vtkIdType i, vtkIdType j, double* k, double t)
134 {
135 if (j < 0)
136 {
137 return 0.0;
138 }
139 if (j == 0)
140 {
141 if (t >= k[i] && t < k[i+1] && k[i] < k[i+1])
142 {
143 return 1.0;
144 }
145 return 0.0;
146 }
147 double term1 = 0.0;
148 if (k[i] < k[i+j])
149 {
150 term1 = (t-k[i])/(k[i+j]-k[i])*N(i, j-1, k, t);
151 }
152 double term2 = 0.0;
153 if (k[i+1] < k[i+j+1])
154 {
155 term2 = (k[i+j+1]-t)/(k[i+j+1]-k[i+1])*N(i+1, j-1, k, t);
156 }
157 return term1 + term2;
158 }
159 #endif
160
161 #if 0
162 double BCubic(vtkIdType i, double* k, double t)
163 {
164 if (t < k[i-2] || t >= k[i+2])
165 {
166 return 0.0;
167 }
168 double temp;
169 if (t >= k[i-2] && t < k[i-1])
170 {
171 temp = 2.0 + t;
172 return temp*temp*temp/6.0;
173 }
174 if (t >= k[i-1] && t < k[i])
175 {
176 return (4 - 6*t*t - 3*t*t*t)/6.0;
177 }
178 if (t >= k[i] && t < k[i+1])
179 {
180 return (4 - 6*t*t + 3*t*t*t)/6.0;
181 }
182 temp = 2.0 - t;
183 return temp*temp*temp;
184 }
185 #endif
186 }
187
vtkSplineGraphEdges()188 vtkSplineGraphEdges::vtkSplineGraphEdges()
189 {
190 this->Spline = vtkCardinalSpline::New();
191 this->XSpline = nullptr;
192 this->YSpline = nullptr;
193 this->ZSpline = nullptr;
194 this->NumberOfSubdivisions = 20;
195 this->SplineType = CUSTOM;
196 }
197
~vtkSplineGraphEdges()198 vtkSplineGraphEdges::~vtkSplineGraphEdges()
199 {
200 if (this->Spline)
201 {
202 this->Spline->Delete();
203 this->Spline = nullptr;
204 }
205 }
206
GetMTime()207 vtkMTimeType vtkSplineGraphEdges::GetMTime()
208 {
209 vtkMTimeType mtime = this->Superclass::GetMTime();
210 if (this->Spline && this->Spline->GetMTime() > mtime)
211 {
212 mtime = this->Spline->GetMTime();
213 }
214 return mtime;
215 }
216
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)217 int vtkSplineGraphEdges::RequestData(
218 vtkInformation *vtkNotUsed(request),
219 vtkInformationVector **inputVector,
220 vtkInformationVector *outputVector)
221 {
222 if (!this->Spline)
223 {
224 vtkErrorMacro("Must have a valid spline.");
225 return 0;
226 }
227
228 // get the info objects
229 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
230 vtkInformation *outInfo = outputVector->GetInformationObject(0);
231
232 // get the input and output
233 vtkGraph *input = vtkGraph::SafeDownCast(
234 inInfo->Get(vtkDataObject::DATA_OBJECT()));
235 vtkGraph *output = vtkGraph::SafeDownCast(
236 outInfo->Get(vtkDataObject::DATA_OBJECT()));
237
238 output->ShallowCopy(input);
239 output->DeepCopyEdgePoints(input);
240
241 if (this->SplineType == CUSTOM)
242 {
243 this->XSpline.TakeReference(this->Spline->NewInstance());
244 this->XSpline->DeepCopy(this->Spline);
245 this->YSpline.TakeReference(this->Spline->NewInstance());
246 this->YSpline->DeepCopy(this->Spline);
247 this->ZSpline.TakeReference(this->Spline->NewInstance());
248 this->ZSpline->DeepCopy(this->Spline);
249 }
250
251 for (vtkIdType i = 0; i < output->GetNumberOfEdges(); ++i)
252 {
253 if (this->SplineType == BSPLINE)
254 {
255 this->GenerateBSpline(output, i);
256 }
257 else
258 {
259 this->GeneratePoints(output, i);
260 }
261 if (i % 1000 == 0)
262 {
263 double progress = static_cast<double>(i)/
264 static_cast<double>(output->GetNumberOfEdges());
265 this->InvokeEvent(vtkCommand::ProgressEvent, &progress);
266 }
267 }
268
269 return 1;
270 }
271
GeneratePoints(vtkGraph * g,vtkIdType e)272 void vtkSplineGraphEdges::GeneratePoints(vtkGraph* g, vtkIdType e)
273 {
274 // Initialize the splines
275 this->XSpline->RemoveAllPoints();
276 this->YSpline->RemoveAllPoints();
277 this->ZSpline->RemoveAllPoints();
278
279 vtkIdType numInternalPoints;
280 double* internalPoints;
281 g->GetEdgePoints(e, numInternalPoints, internalPoints);
282
283 vtkIdType numPoints = numInternalPoints + 2;
284 double* points = new double[3*static_cast<size_t>(numPoints)];
285 memcpy(points + 3, internalPoints, sizeof(double)*3
286 *static_cast<size_t>(numInternalPoints));
287 g->GetPoint(g->GetSourceVertex(e), points);
288 g->GetPoint(g->GetTargetVertex(e), points + 3*(numInternalPoints+1));
289
290 double* xPrev;
291 double* x;
292 double length = 0.0;
293 double* xEnd = points + 3*numPoints;
294 for (xPrev = points, x = points+3; x != xEnd; xPrev += 3, x += 3)
295 {
296 double len = sqrt(vtkMath::Distance2BetweenPoints(x,xPrev));
297 length += len;
298 }
299
300 if (length <= 0.0)
301 {
302 return;
303 }
304
305 // Now we insert points into the splines with the parametric coordinate
306 // based on length. We keep track of the parametric coordinates
307 // of the points for later point interpolation.
308 this->XSpline->AddPoint(0.0, points[0]);
309 this->YSpline->AddPoint(0.0, points[1]);
310 this->ZSpline->AddPoint(0.0, points[2]);
311 double len = 0.0;
312 for (xPrev = points, x = points+3; x != xEnd; xPrev += 3, x += 3)
313 {
314 double dist = sqrt(vtkMath::Distance2BetweenPoints(x, xPrev));
315 if (dist == 0)
316 {
317 continue;
318 }
319 len += dist;
320 double t = len/length;
321
322 this->XSpline->AddPoint(t, x[0]);
323 this->YSpline->AddPoint(t, x[1]);
324 this->ZSpline->AddPoint(t, x[2]);
325 }
326
327 // Now compute the new points
328 vtkIdType numNewPoints = this->NumberOfSubdivisions - 1;
329 double* newPoints = new double[3*numNewPoints];
330 vtkIdType i;
331 for (i = 0, x = newPoints; i < numNewPoints; i++, x += 3)
332 {
333 double t = static_cast<double>(i+1) /
334 static_cast<double>(this->NumberOfSubdivisions);
335 x[0] = this->XSpline->Evaluate(t);
336 x[1] = this->YSpline->Evaluate(t);
337 x[2] = this->ZSpline->Evaluate(t);
338 }
339 g->SetEdgePoints(e, numNewPoints, newPoints);
340
341 delete [] points;
342 delete [] newPoints;
343 }
344
GenerateBSpline(vtkGraph * g,vtkIdType e)345 void vtkSplineGraphEdges::GenerateBSpline(vtkGraph* g, vtkIdType e)
346 {
347 vtkIdType numInternalPoints;
348 double* internalPoints;
349 g->GetEdgePoints(e, numInternalPoints, internalPoints);
350
351 // Duplicate internal point if there is just one, so there are at least
352 // four points, required for B-spline.
353 bool repeat = false;
354 if (numInternalPoints == 1)
355 {
356 repeat = true;
357 numInternalPoints = 2;
358 }
359 vtkIdType numPoints = numInternalPoints + 2;
360 double* points = new double[3*numPoints];
361 if (repeat)
362 {
363 memcpy(points + 3, internalPoints, sizeof(double)*3);
364 memcpy(points + 6, internalPoints, sizeof(double)*3);
365 }
366 else
367 {
368 memcpy(points + 3, internalPoints, sizeof(double)*3
369 *static_cast<size_t>(numInternalPoints));
370 }
371 g->GetPoint(g->GetSourceVertex(e), points);
372 g->GetPoint(g->GetTargetVertex(e), points + 3*(numInternalPoints+1));
373
374 if (numPoints <= 3)
375 {
376 return;
377 }
378
379 // Compute the knot vector
380 vtkIdType numKnots = numPoints + 4;
381 double* knots = new double[numKnots];
382 knots[0] = 0.0;
383 knots[1] = 0.0;
384 knots[2] = 0.0;
385 knots[3] = 0.0;
386 knots[numKnots-4] = 1.0;
387 knots[numKnots-3] = 1.0;
388 knots[numKnots-2] = 1.0;
389 knots[numKnots-1] = 1.0;
390
391 vtkIdType i;
392 for (i = 4; i < numKnots-4; ++i)
393 {
394 knots[i] = static_cast<double>(i-3)/static_cast<double>(numKnots-7);
395 }
396
397 // Special case of 3 points, make symmetric
398 if (numPoints == 3)
399 {
400 knots[3] = 0.5;
401 }
402
403 // Now compute the new points
404 vtkIdType numNewPoints = this->NumberOfSubdivisions - 1;
405 double* newPoints = new double[3*numNewPoints];
406 double* xNew;
407 for (i = 0, xNew = newPoints; i < numNewPoints; i++, xNew += 3)
408 {
409 xNew[0] = 0.0;
410 xNew[1] = 0.0;
411 xNew[2] = 0.0;
412 double t = static_cast<double>(i+1) /
413 static_cast<double>(this->NumberOfSubdivisions);
414 double* x = points;
415 //double bsum = 0.0;
416 for (vtkIdType j = 0; j < numPoints; ++j, x += 3)
417 {
418 //double b = BCubic(j+2, knots, t);
419 //double b = N(j, 3, knots, t);
420 double b = CubicSpline(j, knots, t);
421 //bsum += b;
422 xNew[0] += x[0]*b;
423 xNew[1] += x[1]*b;
424 xNew[2] += x[2]*b;
425 }
426 //cerr << "bsum: " << bsum << endl;
427 }
428 g->SetEdgePoints(e, numNewPoints, newPoints);
429
430 delete [] points;
431 delete [] knots;
432 delete [] newPoints;
433 }
434
PrintSelf(ostream & os,vtkIndent indent)435 void vtkSplineGraphEdges::PrintSelf(ostream& os, vtkIndent indent)
436 {
437 this->Superclass::PrintSelf(os,indent);
438 os << indent << "SplineType: " << this->SplineType << endl;
439 os << indent << "NumberOfSubdivisions: " << this->NumberOfSubdivisions << endl;
440 os << indent << "Spline: " << (this->Spline ? "" : "(none)") << endl;
441 if (this->Spline)
442 {
443 this->Spline->PrintSelf(os, indent.GetNextIndent());
444 }
445 }
446