1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkParametricFunctionSource.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 "vtkParametricFunctionSource.h"
16 #include "vtkCellArray.h"
17 #include "vtkFloatArray.h"
18 #include "vtkInformation.h"
19 #include "vtkInformationVector.h"
20 #include "vtkMath.h"
21 #include "vtkObjectFactory.h"
22 #include "vtkParametricFunction.h"
23 #include "vtkPointData.h"
24 #include "vtkPoints.h"
25 #include "vtkPolyData.h"
26 #include "vtkPolyDataNormals.h"
27 #include "vtkSmartPointer.h"
28 #include "vtkTriangleFilter.h"
29
30 #include <cmath>
31 #include <string>
32
33 vtkStandardNewMacro(vtkParametricFunctionSource);
34 vtkCxxSetObjectMacro(vtkParametricFunctionSource, ParametricFunction, vtkParametricFunction);
35
36 //------------------------------------------------------------------------------
vtkParametricFunctionSource()37 vtkParametricFunctionSource::vtkParametricFunctionSource()
38 : ParametricFunction(nullptr)
39 , UResolution(50)
40 , VResolution(50)
41 , WResolution(50)
42 , GenerateTextureCoordinates(0)
43 , ScalarMode(vtkParametricFunctionSource::SCALAR_NONE)
44 , OutputPointsPrecision(vtkAlgorithm::SINGLE_PRECISION)
45 {
46 this->SetNumberOfInputPorts(0);
47 this->GenerateNormals = 1;
48 }
49
50 //------------------------------------------------------------------------------
~vtkParametricFunctionSource()51 vtkParametricFunctionSource::~vtkParametricFunctionSource()
52 {
53 this->SetParametricFunction(nullptr);
54 }
55
56 namespace
57 {
58 /**
59 * Make the cells containing the ordered point Ids.
60 *
61 */
AddTriCells(vtkCellArray * cellArray,int id1,int id2,int id3,int id4,bool clockwise)62 void AddTriCells(vtkCellArray* cellArray, int id1, int id2, int id3, int id4, bool clockwise)
63 {
64 cellArray->InsertNextCell(3);
65 if (clockwise)
66 {
67 cellArray->InsertCellPoint(id1);
68 cellArray->InsertCellPoint(id2);
69 cellArray->InsertCellPoint(id3);
70 cellArray->InsertNextCell(3);
71 cellArray->InsertCellPoint(id1);
72 cellArray->InsertCellPoint(id3);
73 cellArray->InsertCellPoint(id4);
74 }
75 else
76 {
77 cellArray->InsertCellPoint(id1);
78 cellArray->InsertCellPoint(id3);
79 cellArray->InsertCellPoint(id2);
80 cellArray->InsertNextCell(3);
81 cellArray->InsertCellPoint(id1);
82 cellArray->InsertCellPoint(id4);
83 cellArray->InsertCellPoint(id3);
84 }
85 }
86
87 } // anonymous namespace
88
89 //------------------------------------------------------------------------------
MakeTriangles(vtkCellArray * cells,int PtsU,int PtsV)90 void vtkParametricFunctionSource::MakeTriangles(vtkCellArray* cells, int PtsU, int PtsV)
91 {
92 int id1 = 0;
93 int id2 = 0;
94 int id3 = 0;
95 int id4 = 0;
96
97 vtkDebugMacro(<< "Executing MakeTriangles()");
98
99 bool clockwise = (this->ParametricFunction->GetClockwiseOrdering() != 0);
100
101 vtkIdType numCells = (PtsU + this->ParametricFunction->GetJoinU() - 1) *
102 (PtsV + this->ParametricFunction->GetJoinV() - 1) * 2;
103 cells->AllocateExact(numCells, numCells * 3);
104
105 for (int i = 0; i < PtsU - 1; ++i)
106 {
107 // Fill the allocated space with the indexes to the points.
108 for (int j = 0; j < PtsV - 1; ++j)
109 {
110 id1 = j + i * PtsV;
111 id2 = id1 + PtsV;
112 id3 = id2 + 1;
113 id4 = id1 + 1;
114 AddTriCells(cells, id1, id2, id3, id4, clockwise);
115 }
116 // If necessary, connect the ends of the triangle strip.
117 if (this->ParametricFunction->GetJoinV())
118 {
119 id1 = id4;
120 id2 = id3;
121 if (this->ParametricFunction->GetTwistV())
122 {
123 id3 = (i + 1) * PtsV;
124 id4 = i * PtsV;
125 }
126 else
127 {
128 id3 = i * PtsV;
129 id4 = (i + 1) * PtsV;
130 }
131 AddTriCells(cells, id1, id2, id3, id4, clockwise);
132 }
133 }
134 // If required, connect the last triangle strip to the first by
135 // adding a new triangle strip and filling it with the indexes
136 // to the points.
137 if (this->ParametricFunction->GetJoinU())
138 {
139 for (int j = 0; j < PtsV - 1; ++j)
140 {
141 id1 = j + (PtsU - 1) * PtsV;
142 id3 = id1 + 1;
143 if (this->ParametricFunction->GetTwistU())
144 {
145 id2 = PtsV - 1 - j;
146 id4 = id2 - 1;
147 }
148 else
149 {
150 id2 = j;
151 id4 = id2 + 1;
152 }
153 AddTriCells(cells, id1, id2, id3, id4, clockwise);
154 }
155
156 // If necessary, connect the ends of the triangle strip.
157 if (this->ParametricFunction->GetJoinV())
158 {
159 id1 = id3;
160 id2 = id4;
161 if (this->ParametricFunction->GetTwistU())
162 {
163 if (this->ParametricFunction->GetTwistV())
164 {
165 id3 = PtsV - 1;
166 id4 = (PtsU - 1) * PtsV;
167 }
168 else
169 {
170 id3 = (PtsU - 1) * PtsV;
171 id4 = PtsV - 1;
172 }
173 }
174 else
175 {
176 if (this->ParametricFunction->GetTwistV())
177 {
178 id3 = 0;
179 id4 = (PtsU - 1) * PtsV;
180 }
181 else
182 {
183 id3 = (PtsU - 1) * PtsV;
184 id4 = 0;
185 }
186 }
187 AddTriCells(cells, id1, id2, id3, id4, clockwise);
188 }
189 }
190 cells->Modified();
191 vtkDebugMacro(<< "MakeTriangles() finished.");
192 }
193
194 //------------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (info),vtkInformationVector ** vtkNotUsed (inputV),vtkInformationVector * output)195 int vtkParametricFunctionSource::RequestData(vtkInformation* vtkNotUsed(info),
196 vtkInformationVector** vtkNotUsed(inputV), vtkInformationVector* output)
197 {
198 vtkDebugMacro(<< "Executing");
199
200 // Check that a parametric function has been defined
201 if (!this->ParametricFunction)
202 {
203 vtkErrorMacro(<< "Parametric function not defined");
204 return 1;
205 }
206
207 switch (this->ParametricFunction->GetDimension())
208 {
209 case 1:
210 this->Produce1DOutput(output);
211 break;
212 case 2:
213 this->Produce2DOutput(output);
214 break;
215 default:
216 vtkErrorMacro("Functions of dimension " << this->ParametricFunction->GetDimension()
217 << " are not supported.");
218 }
219
220 return 1;
221 }
222
223 //------------------------------------------------------------------------------
Produce1DOutput(vtkInformationVector * output)224 void vtkParametricFunctionSource::Produce1DOutput(vtkInformationVector* output)
225 {
226 vtkIdType numPts = this->UResolution + 1;
227 vtkSmartPointer<vtkCellArray> lines = vtkSmartPointer<vtkCellArray>::New();
228 vtkSmartPointer<vtkPoints> pts = vtkSmartPointer<vtkPoints>::New();
229
230 // Set the desired precision for the points in the output.
231 if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
232 {
233 pts->SetDataType(VTK_DOUBLE);
234 }
235 else
236 {
237 pts->SetDataType(VTK_FLOAT);
238 }
239
240 pts->SetNumberOfPoints(numPts);
241 vtkIdType i;
242 double x[3], Du[3], t[3];
243
244 lines->AllocateEstimate(1, numPts);
245 lines->InsertNextCell(numPts);
246
247 // Insert points and cell points
248 for (i = 0; i < numPts; i++)
249 {
250 t[0] = (double)i / this->UResolution;
251 this->ParametricFunction->Evaluate(t, x, Du);
252 pts->SetPoint(i, x);
253 lines->InsertCellPoint(i);
254 }
255
256 vtkInformation* outInfo = output->GetInformationObject(0);
257 vtkPolyData* outData = static_cast<vtkPolyData*>(outInfo->Get(vtkDataObject::DATA_OBJECT()));
258 outData->SetPoints(pts);
259 outData->SetLines(lines);
260 }
261
262 //------------------------------------------------------------------------------
Produce2DOutput(vtkInformationVector * output)263 void vtkParametricFunctionSource::Produce2DOutput(vtkInformationVector* output)
264 {
265 // Adjust so the ranges:
266 // this->MinimumU ... this->ParametricFunction->GetMaximumU(),
267 // this->MinimumV ... this->ParametricFunction->GetMaximumV()
268 // are included in the triangulation.
269 double MaxU = this->ParametricFunction->GetMaximumU() +
270 (this->ParametricFunction->GetMaximumU() - this->ParametricFunction->GetMinimumU()) /
271 (this->UResolution - 1);
272 int PtsU = this->UResolution;
273 double MaxV = this->ParametricFunction->GetMaximumV() +
274 (this->ParametricFunction->GetMaximumV() - this->ParametricFunction->GetMinimumV()) /
275 (this->VResolution - 1);
276 int PtsV = this->VResolution;
277 int totPts = PtsU * PtsV;
278
279 // Scalars associated with each point
280 vtkSmartPointer<vtkFloatArray> sval = vtkSmartPointer<vtkFloatArray>::New();
281 if (this->ScalarMode != SCALAR_NONE)
282 {
283 sval->SetNumberOfTuples(totPts);
284 sval->SetName("Scalars");
285 }
286
287 // The normals to the surface
288 vtkSmartPointer<vtkFloatArray> nval = vtkSmartPointer<vtkFloatArray>::New();
289 if (this->GenerateNormals)
290 {
291 nval->SetNumberOfComponents(3);
292 nval->SetNumberOfTuples(totPts);
293 nval->SetName("Normals");
294 }
295
296 // Texture coordinates
297 vtkSmartPointer<vtkFloatArray> newTCoords = vtkSmartPointer<vtkFloatArray>::New();
298 if (this->GenerateTextureCoordinates != 0)
299 {
300 newTCoords->SetNumberOfComponents(2);
301 newTCoords->Allocate(2 * totPts);
302 newTCoords->SetName("Textures");
303 }
304
305 vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
306
307 // Set the desired precision for the points in the output.
308 if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
309 {
310 points->SetDataType(VTK_DOUBLE);
311 }
312 else
313 {
314 points->SetDataType(VTK_FLOAT);
315 }
316
317 points->SetNumberOfPoints(totPts);
318
319 double uStep = (MaxU - this->ParametricFunction->GetMinimumU()) / PtsU;
320 double vStep = (MaxV - this->ParametricFunction->GetMinimumV()) / PtsV;
321
322 // Find the mid points of the (u,v) map.
323 double u0 = this->ParametricFunction->GetMinimumU();
324 double u_mp = (MaxU - u0) / 2.0 + u0 - uStep;
325 while (u0 < u_mp)
326 {
327 u0 += uStep;
328 }
329
330 double v0 = this->ParametricFunction->GetMinimumV();
331 double v_mp = (MaxV - v0) / 2.0 + v0 - vStep;
332 while (v0 < v_mp)
333 {
334 v0 += vStep;
335 }
336 u_mp += uStep;
337 v_mp += vStep;
338
339 // At this point (u_mp, v_mp) is the midpoint of the (u,v) map and (u0,v0)
340 // corresponds to the nearest grid point to the midpoint of the (u,v) map.
341 //
342 double rel_u = 0; // will be u - u_mp
343 double rel_v = 0; // will be v - v_mp
344
345 int k = 0;
346 double uv[3];
347 uv[0] = this->ParametricFunction->GetMinimumU() - uStep;
348
349 float MaxI = PtsU - 1;
350 float MaxJ = PtsV - 1;
351
352 for (int i = 0; i < PtsU; ++i)
353 {
354 uv[0] += uStep;
355 uv[1] = this->ParametricFunction->GetMinimumV() - vStep;
356
357 double tc[2];
358 if (this->GenerateTextureCoordinates != 0)
359 {
360 tc[0] = i / MaxI;
361 }
362
363 for (int j = 0; j < PtsV; ++j)
364 {
365 uv[1] += vStep;
366
367 if (this->GenerateTextureCoordinates != 0)
368 {
369 tc[1] = 1.0 - j / MaxJ;
370 newTCoords->InsertNextTuple(tc);
371 }
372
373 // The point
374 double Pt[3];
375 // Partial derivative at Pt with respect to u,v,w.
376 double Du[9];
377 // Partial derivative at Pt with respect to v.
378 double* Dv = Du + 3;
379
380 // Calculate fn(u)->(Pt,Du).
381 this->ParametricFunction->Evaluate(uv, Pt, Du);
382
383 // Insert the points and scalar.
384 points->InsertPoint(k, Pt[0], Pt[1], Pt[2]);
385 double scalar;
386
387 if (this->ScalarMode != SCALAR_NONE)
388 {
389 switch (this->ScalarMode)
390 {
391 case SCALAR_U:
392 scalar = uv[0];
393 break;
394 case SCALAR_V:
395 scalar = uv[1];
396 break;
397 case SCALAR_U0:
398 scalar = uv[0] == u0 ? 1 : 0;
399 break;
400 case SCALAR_V0:
401 scalar = uv[1] == v0 ? 1 : 0;
402 break;
403 case SCALAR_U0V0:
404 scalar = 0;
405 // u0, v0
406 if (uv[0] == u0 && uv[1] == v0)
407 {
408 scalar = 3;
409 }
410 else
411 {
412 // u0 line
413 if (uv[0] == u0)
414 {
415 scalar = 1;
416 }
417 else
418 {
419 // v0 line
420 if (uv[1] == v0)
421 {
422 scalar = 2;
423 }
424 }
425 }
426 break;
427 case SCALAR_MODULUS:
428 rel_u = uv[0] - u_mp;
429 rel_v = uv[1] - v_mp;
430 scalar = sqrt(rel_u * rel_u + rel_v * rel_v);
431 break;
432 case SCALAR_PHASE:
433 rel_u = uv[0] - u_mp;
434 rel_v = uv[1] - v_mp;
435 if (rel_v == 0 && rel_u == 0)
436 {
437 scalar = 0;
438 }
439 else
440 {
441 scalar = vtkMath::DegreesFromRadians(atan2(rel_v, rel_u));
442 if (scalar < 0)
443 {
444 scalar += 360;
445 }
446 }
447 break;
448 case SCALAR_QUADRANT:
449 if (uv[0] >= u0 && uv[1] >= v0)
450 {
451 scalar = 1;
452 break;
453 }
454 if (uv[0] < u0 && uv[1] >= v0)
455 {
456 scalar = 2;
457 break;
458 }
459 if (uv[0] < u0 && uv[1] < v0)
460 {
461 scalar = 3;
462 }
463 else
464 {
465 scalar = 4;
466 }
467 break;
468 case SCALAR_X:
469 scalar = Pt[0];
470 break;
471 case SCALAR_Y:
472 scalar = Pt[1];
473 break;
474 case SCALAR_Z:
475 scalar = Pt[2];
476 break;
477 case SCALAR_DISTANCE:
478 scalar = sqrt(Pt[0] * Pt[0] + Pt[1] * Pt[1] + Pt[2] * Pt[2]);
479 break;
480 case SCALAR_FUNCTION_DEFINED:
481 scalar = this->ParametricFunction->EvaluateScalar(uv, Pt, Du);
482 break;
483 case SCALAR_NONE:
484 default:
485 scalar = 0;
486 }
487 sval->SetValue(k, scalar);
488 }
489
490 // Calculate the normal.
491 if (this->ParametricFunction->GetDerivativesAvailable() && this->GenerateNormals)
492 {
493 double n[3];
494 if (this->ParametricFunction->GetClockwiseOrdering() == 0)
495 {
496 // Anti-clockwise ordering
497 vtkMath::Cross(Dv, Du, n);
498 }
499 else
500 {
501 // Clockwise ordering
502 vtkMath::Cross(Du, Dv, n);
503 }
504 nval->SetTuple3(k, n[0], n[1], n[2]);
505 }
506
507 ++k;
508 }
509 }
510
511 vtkInformation* outInfo = output->GetInformationObject(0);
512 vtkPolyData* outData = static_cast<vtkPolyData*>(outInfo->Get(vtkDataObject::DATA_OBJECT()));
513 vtkCellArray* tris = vtkCellArray::New();
514 this->MakeTriangles(tris, PtsU, PtsV);
515 outData->SetPoints(points);
516 outData->SetPolys(tris);
517
518 if (this->GenerateNormals)
519 {
520 if (this->ParametricFunction->GetDerivativesAvailable())
521 {
522 outData->GetPointData()->SetNormals(nval);
523 }
524 else
525 {
526 // Used to hold the surface
527 vtkSmartPointer<vtkPolyData> pd = vtkSmartPointer<vtkPolyData>::New();
528 pd->SetPoints(points);
529 pd->SetPolys(tris);
530 vtkSmartPointer<vtkPolyDataNormals> norm = vtkSmartPointer<vtkPolyDataNormals>::New();
531 // we prevent vtkPolyDataNormals from generating new points
532 // so that the number of newTCoords matches the number of points.
533 norm->SplittingOff();
534 norm->SetInputData(pd);
535 norm->Update();
536 outData->DeepCopy(norm->GetOutput());
537 }
538 }
539
540 tris->Delete();
541 if (this->ScalarMode != SCALAR_NONE)
542 {
543 outData->GetPointData()->SetScalars(sval);
544 }
545 if (this->GenerateTextureCoordinates != 0)
546 {
547 outData->GetPointData()->SetTCoords(newTCoords);
548 }
549 outData->Modified();
550 }
551
552 //------------------------------------------------------------------------------
GetMTime()553 vtkMTimeType vtkParametricFunctionSource::GetMTime()
554 {
555 vtkMTimeType mTime = this->Superclass::GetMTime();
556 vtkMTimeType funcMTime;
557
558 if (this->ParametricFunction != nullptr)
559 {
560 funcMTime = this->ParametricFunction->GetMTime();
561 mTime = (funcMTime > mTime ? funcMTime : mTime);
562 }
563
564 return mTime;
565 }
566
567 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)568 void vtkParametricFunctionSource::PrintSelf(ostream& os, vtkIndent indent)
569 {
570 this->Superclass::PrintSelf(os, indent);
571
572 os << indent << "U Resolution: " << this->UResolution << "\n";
573 os << indent << "V Resolution: " << this->VResolution << "\n";
574 os << indent << "W Resolution: " << this->WResolution << "\n";
575
576 if (this->ParametricFunction)
577 {
578 os << indent << "Parametric Function: " << this->ParametricFunction << "\n";
579 }
580 else
581 {
582 os << indent << "No Parametric function defined\n";
583 }
584
585 std::string s;
586 switch (this->ScalarMode)
587 {
588 case SCALAR_NONE:
589 s = "SCALAR_NONE";
590 break;
591 case SCALAR_U:
592 s = "SCALAR_U";
593 break;
594 case SCALAR_V:
595 s = "SCALAR_V";
596 break;
597 case SCALAR_U0:
598 s = "SCALAR_U0";
599 break;
600 case SCALAR_V0:
601 s = "SCALAR_V0";
602 break;
603 case SCALAR_U0V0:
604 s = "SCALAR_U0V0";
605 break;
606 case SCALAR_MODULUS:
607 s = "SCALAR_MODULUS";
608 break;
609 case SCALAR_PHASE:
610 s = "SCALAR_PHASE";
611 break;
612 case SCALAR_QUADRANT:
613 s = "SCALAR_QUADRANT";
614 break;
615 case SCALAR_X:
616 s = "SCALAR_X";
617 break;
618 case SCALAR_Y:
619 s = "SCALAR_Y";
620 break;
621 case SCALAR_Z:
622 s = "SCALAR_Z";
623 break;
624 case SCALAR_DISTANCE:
625 s = "SCALAR_DISTANCE";
626 break;
627 case SCALAR_FUNCTION_DEFINED:
628 s = "SCALAR_FUNCTION_DEFINED";
629 break;
630 default:
631 s = "Unknown scalar mode.";
632 }
633 os << indent << "Scalar Mode: " << s.c_str() << "\n";
634 os << indent << "GenerateTextureCoordinates:" << (this->GenerateTextureCoordinates ? "On" : "Off")
635 << "\n";
636 os << indent << "Output Points Precision: " << this->OutputPointsPrecision << "\n";
637 }
638