1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkTensorGlyph.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 "vtkTensorGlyph.h"
16
17
18 #include "vtkStreamingDemandDrivenPipeline.h"
19 #include "vtkCell.h"
20 #include "vtkCellArray.h"
21 #include "vtkDataSet.h"
22 #include "vtkExecutive.h"
23 #include "vtkFloatArray.h"
24 #include "vtkMath.h"
25 #include "vtkInformation.h"
26 #include "vtkInformationVector.h"
27 #include "vtkObjectFactory.h"
28 #include "vtkPointData.h"
29 #include "vtkPolyData.h"
30 #include "vtkTransform.h"
31
32 vtkStandardNewMacro(vtkTensorGlyph);
33
34 // Construct object with scaling on and scale factor 1.0. Eigenvalues are
35 // extracted, glyphs are colored with input scalar data, and logarithmic
36 // scaling is turned off.
vtkTensorGlyph()37 vtkTensorGlyph::vtkTensorGlyph()
38 {
39 this->Scaling = 1;
40 this->ScaleFactor = 1.0;
41 this->ExtractEigenvalues = 1;
42 this->ColorGlyphs = 1;
43 this->ColorMode = COLOR_BY_SCALARS;
44 this->ClampScaling = 0;
45 this->MaxScaleFactor = 100;
46 this->ThreeGlyphs = 0;
47 this->Symmetric = 0;
48 this->Length = 1.0;
49
50 this->SetNumberOfInputPorts(2);
51
52 // by default, process active point tensors
53 this->SetInputArrayToProcess(0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS,
54 vtkDataSetAttributes::TENSORS);
55
56 // by default, process active point scalars
57 this->SetInputArrayToProcess(1, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS,
58 vtkDataSetAttributes::SCALARS);
59 }
60
61 //----------------------------------------------------------------------------
62 vtkTensorGlyph::~vtkTensorGlyph() = default;
63
64 //----------------------------------------------------------------------------
RequestUpdateExtent(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)65 int vtkTensorGlyph::RequestUpdateExtent(
66 vtkInformation *vtkNotUsed(request),
67 vtkInformationVector **inputVector,
68 vtkInformationVector *outputVector)
69 {
70 // get the info objects
71 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
72 vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
73 vtkInformation *outInfo = outputVector->GetInformationObject(0);
74
75 if (sourceInfo)
76 {
77 sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
78 0);
79 sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
80 1);
81 sourceInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
82 0);
83 }
84
85 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER(),
86 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_PIECE_NUMBER()));
87 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES(),
88 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_PIECES()));
89 inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
90 outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()));
91 inInfo->Set(vtkStreamingDemandDrivenPipeline::EXACT_EXTENT(), 1);
92
93 return 1;
94 }
95
96 //----------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)97 int vtkTensorGlyph::RequestData(
98 vtkInformation *vtkNotUsed(request),
99 vtkInformationVector **inputVector,
100 vtkInformationVector *outputVector)
101 {
102 // get the info objects
103 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
104 vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
105 vtkInformation *outInfo = outputVector->GetInformationObject(0);
106
107 // get the input and output
108 vtkDataSet *input = vtkDataSet::SafeDownCast(
109 inInfo->Get(vtkDataObject::DATA_OBJECT()));
110 vtkPolyData *source = vtkPolyData::SafeDownCast(
111 sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
112 vtkPolyData *output = vtkPolyData::SafeDownCast(
113 outInfo->Get(vtkDataObject::DATA_OBJECT()));
114
115 vtkDataArray *inTensors;
116 double tensor[9];
117 vtkDataArray *inScalars;
118 vtkIdType numPts, numSourcePts, numSourceCells, inPtId, i;
119 int j;
120 vtkPoints *sourcePts;
121 vtkDataArray *sourceNormals;
122 vtkCellArray *sourceCells, *cells;
123 vtkPoints *newPts;
124 vtkFloatArray *newScalars=nullptr;
125 vtkFloatArray *newNormals=nullptr;
126 double x[3], s;
127 vtkTransform *trans;
128 vtkCell *cell;
129 vtkIdList *cellPts;
130 int npts;
131 vtkIdType *pts;
132 vtkIdType ptIncr, cellId;
133 vtkIdType subIncr;
134 int numDirs, dir, eigen_dir, symmetric_dir;
135 vtkMatrix4x4 *matrix;
136 double *m[3], w[3], *v[3];
137 double m0[3], m1[3], m2[3];
138 double v0[3], v1[3], v2[3];
139 double xv[3], yv[3], zv[3];
140 double maxScale;
141
142 numDirs = (this->ThreeGlyphs?3:1)*(this->Symmetric+1);
143
144 // set up working matrices
145 m[0] = m0; m[1] = m1; m[2] = m2;
146 v[0] = v0; v[1] = v1; v[2] = v2;
147
148 vtkDebugMacro(<<"Generating tensor glyphs");
149
150 vtkPointData *outPD = output->GetPointData();
151 inTensors = this->GetInputArrayToProcess(0, inputVector);
152 inScalars = this->GetInputArrayToProcess(1, inputVector);
153 numPts = input->GetNumberOfPoints();
154
155 if ( !inTensors || numPts < 1 )
156 {
157 vtkErrorMacro(<<"No data to glyph!");
158 return 1;
159 }
160
161 pts = new vtkIdType[source->GetMaxCellSize()];
162 trans = vtkTransform::New();
163 matrix = vtkMatrix4x4::New();
164
165 //
166 // Allocate storage for output PolyData
167 //
168 sourcePts = source->GetPoints();
169 numSourcePts = sourcePts->GetNumberOfPoints();
170 numSourceCells = source->GetNumberOfCells();
171
172 newPts = vtkPoints::New();
173 newPts->Allocate(numDirs*numPts*numSourcePts);
174
175 // Setting up for calls to PolyData::InsertNextCell()
176 if ( (sourceCells=source->GetVerts())->GetNumberOfCells() > 0 )
177 {
178 cells = vtkCellArray::New();
179 cells->Allocate(numDirs*numPts*sourceCells->GetSize());
180 output->SetVerts(cells);
181 cells->Delete();
182 }
183 if ( (sourceCells=this->GetSource()->GetLines())->GetNumberOfCells() > 0 )
184 {
185 cells = vtkCellArray::New();
186 cells->Allocate(numDirs*numPts*sourceCells->GetSize());
187 output->SetLines(cells);
188 cells->Delete();
189 }
190 if ( (sourceCells=this->GetSource()->GetPolys())->GetNumberOfCells() > 0 )
191 {
192 cells = vtkCellArray::New();
193 cells->Allocate(numDirs*numPts*sourceCells->GetSize());
194 output->SetPolys(cells);
195 cells->Delete();
196 }
197 if ( (sourceCells=this->GetSource()->GetStrips())->GetNumberOfCells() > 0 )
198 {
199 cells = vtkCellArray::New();
200 cells->Allocate(numDirs*numPts*sourceCells->GetSize());
201 output->SetStrips(cells);
202 cells->Delete();
203 }
204
205 // only copy scalar data through
206 vtkPointData *pd = this->GetSource()->GetPointData();
207 // generate scalars if eigenvalues are chosen or if scalars exist.
208 if (this->ColorGlyphs &&
209 ((this->ColorMode == COLOR_BY_EIGENVALUES) ||
210 (inScalars && (this->ColorMode == COLOR_BY_SCALARS)) ) )
211 {
212 newScalars = vtkFloatArray::New();
213 newScalars->Allocate(numDirs*numPts*numSourcePts);
214 if (this->ColorMode == COLOR_BY_EIGENVALUES)
215 {
216 newScalars->SetName("MaxEigenvalue");
217 }
218 else
219 {
220 newScalars->SetName(inScalars->GetName());
221 }
222 }
223 else
224 {
225 outPD->CopyAllOff();
226 outPD->CopyScalarsOn();
227 outPD->CopyAllocate(pd,numDirs*numPts*numSourcePts);
228 }
229 if ( (sourceNormals = pd->GetNormals()) )
230 {
231 newNormals = vtkFloatArray::New();
232 newNormals->SetNumberOfComponents(3);
233 newNormals->SetName("Normals");
234 newNormals->Allocate(numDirs*3*numPts*numSourcePts);
235 }
236 //
237 // First copy all topology (transformation independent)
238 //
239 for (inPtId=0; inPtId < numPts; inPtId++)
240 {
241 ptIncr = numDirs * inPtId * numSourcePts;
242 for (cellId=0; cellId < numSourceCells; cellId++)
243 {
244 cell = this->GetSource()->GetCell(cellId);
245 cellPts = cell->GetPointIds();
246 npts = cellPts->GetNumberOfIds();
247 for (dir=0; dir < numDirs; dir++)
248 {
249 // This variable may be removed, but that
250 // will not improve readability
251 subIncr = ptIncr + dir*numSourcePts;
252 for (i=0; i < npts; i++)
253 {
254 pts[i] = cellPts->GetId(i) + subIncr;
255 }
256 output->InsertNextCell(cell->GetCellType(),npts,pts);
257 }
258 }
259 }
260 //
261 // Traverse all Input points, transforming glyph at Source points
262 //
263 trans->PreMultiply();
264
265 for (inPtId=0; inPtId < numPts; inPtId++)
266 {
267 ptIncr = numDirs * inPtId * numSourcePts;
268
269 // Translation is postponed
270 // Symmetric tensor support
271 inTensors->GetTuple(inPtId, tensor);
272 if (inTensors->GetNumberOfComponents() == 6)
273 {
274 vtkMath::TensorFromSymmetricTensor(tensor);
275 }
276
277 // compute orientation vectors and scale factors from tensor
278 if ( this->ExtractEigenvalues ) // extract appropriate eigenfunctions
279 {
280 // We are interested in the symmetrical part of the tensor only, since
281 // eigenvalues are real if and only if the matrice of reals is symmetrical
282 for (j=0; j<3; j++)
283 {
284 for (i=0; i<3; i++)
285 {
286 m[i][j] = 0.5 * (tensor[i + 3 * j] + tensor[j + 3 * i]);
287 }
288 }
289 vtkMath::Jacobi(m, w, v);
290
291 //copy eigenvectors
292 xv[0] = v[0][0]; xv[1] = v[1][0]; xv[2] = v[2][0];
293 yv[0] = v[0][1]; yv[1] = v[1][1]; yv[2] = v[2][1];
294 zv[0] = v[0][2]; zv[1] = v[1][2]; zv[2] = v[2][2];
295 }
296 else //use tensor columns as eigenvectors
297 {
298 for (i=0; i<3; i++)
299 {
300 xv[i] = tensor[i];
301 yv[i] = tensor[i+3];
302 zv[i] = tensor[i+6];
303 }
304 w[0] = vtkMath::Normalize(xv);
305 w[1] = vtkMath::Normalize(yv);
306 w[2] = vtkMath::Normalize(zv);
307 }
308
309 // compute scale factors
310 w[0] *= this->ScaleFactor;
311 w[1] *= this->ScaleFactor;
312 w[2] *= this->ScaleFactor;
313
314 if ( this->ClampScaling )
315 {
316 for (maxScale=0.0, i=0; i<3; i++)
317 {
318 if ( maxScale < fabs(w[i]) )
319 {
320 maxScale = fabs(w[i]);
321 }
322 }
323 if ( maxScale > this->MaxScaleFactor )
324 {
325 maxScale = this->MaxScaleFactor / maxScale;
326 for (i=0; i<3; i++)
327 {
328 w[i] *= maxScale; //preserve overall shape of glyph
329 }
330 }
331 }
332
333 // normalization is postponed
334
335 // make sure scale is okay (non-zero) and scale data
336 for (maxScale=0.0, i=0; i<3; i++)
337 {
338 if ( w[i] > maxScale )
339 {
340 maxScale = w[i];
341 }
342 }
343 if ( maxScale == 0.0 )
344 {
345 maxScale = 1.0;
346 }
347 for (i=0; i<3; i++)
348 {
349 if ( w[i] == 0.0 )
350 {
351 w[i] = maxScale * 1.0e-06;
352 }
353 }
354
355 // Now do the real work for each "direction"
356
357 for (dir=0; dir < numDirs; dir++)
358 {
359 eigen_dir = dir%(this->ThreeGlyphs?3:1);
360 symmetric_dir = dir/(this->ThreeGlyphs?3:1);
361
362 // Remove previous scales ...
363 trans->Identity();
364
365 // translate Source to Input point
366 input->GetPoint(inPtId, x);
367 trans->Translate(x[0], x[1], x[2]);
368
369 // normalized eigenvectors rotate object for eigen direction 0
370 matrix->Element[0][0] = xv[0];
371 matrix->Element[0][1] = yv[0];
372 matrix->Element[0][2] = zv[0];
373 matrix->Element[1][0] = xv[1];
374 matrix->Element[1][1] = yv[1];
375 matrix->Element[1][2] = zv[1];
376 matrix->Element[2][0] = xv[2];
377 matrix->Element[2][1] = yv[2];
378 matrix->Element[2][2] = zv[2];
379 trans->Concatenate(matrix);
380
381 if (eigen_dir == 1)
382 {
383 trans->RotateZ(90.0);
384 }
385
386 if (eigen_dir == 2)
387 {
388 trans->RotateY(-90.0);
389 }
390
391 if (this->ThreeGlyphs)
392 {
393 trans->Scale(w[eigen_dir], this->ScaleFactor, this->ScaleFactor);
394 }
395 else
396 {
397 trans->Scale(w[0], w[1], w[2]);
398 }
399
400 // Mirror second set to the symmetric position
401 if (symmetric_dir == 1)
402 {
403 trans->Scale(-1.,1.,1.);
404 }
405
406 // if the eigenvalue is negative, shift to reverse direction.
407 // The && is there to ensure that we do not change the
408 // old behaviour of vtkTensorGlyphs (which only used one dir),
409 // in case there is an oriented glyph, e.g. an arrow.
410 if (w[eigen_dir] < 0 && numDirs > 1)
411 {
412 trans->Translate(-this->Length, 0., 0.);
413 }
414
415 // multiply points (and normals if available) by resulting
416 // matrix
417 trans->TransformPoints(sourcePts,newPts);
418
419 // Apply the transformation to a series of points,
420 // and append the results to outPts.
421 if ( newNormals )
422 {
423 // a negative determinant means the transform turns the
424 // glyph surface inside out, and its surface normals all
425 // point inward. The following scale corrects the surface
426 // normals to point outward.
427 if (trans->GetMatrix()->Determinant() < 0)
428 {
429 trans->Scale(-1.0,-1.0,-1.0);
430 }
431 trans->TransformNormals(sourceNormals,newNormals);
432 }
433
434 // Copy point data from source
435 if ( this->ColorGlyphs && inScalars &&
436 (this->ColorMode == COLOR_BY_SCALARS) )
437 {
438 s = inScalars->GetComponent(inPtId, 0);
439 for (i=0; i < numSourcePts; i++)
440 {
441 newScalars->InsertTuple(ptIncr+i, &s);
442 }
443 }
444 else if (this->ColorGlyphs &&
445 (this->ColorMode == COLOR_BY_EIGENVALUES) )
446 {
447 // If ThreeGlyphs is false we use the first (largest)
448 // eigenvalue as scalar.
449 s = w[eigen_dir];
450 for (i=0; i < numSourcePts; i++)
451 {
452 newScalars->InsertTuple(ptIncr+i, &s);
453 }
454 }
455 else
456 {
457 for (i=0; i < numSourcePts; i++)
458 {
459 outPD->CopyData(pd,i,ptIncr+i);
460 }
461 }
462 ptIncr += numSourcePts;
463 }
464 }
465 vtkDebugMacro(<<"Generated " << numPts <<" tensor glyphs");
466 //
467 // Update output and release memory
468 //
469 delete [] pts;
470
471 output->SetPoints(newPts);
472 newPts->Delete();
473
474 if ( newScalars )
475 {
476 int idx = outPD->AddArray(newScalars);
477 outPD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
478 newScalars->Delete();
479 }
480
481 if ( newNormals )
482 {
483 outPD->SetNormals(newNormals);
484 newNormals->Delete();
485 }
486
487 output->Squeeze();
488 trans->Delete();
489 matrix->Delete();
490
491 return 1;
492 }
493
494 //----------------------------------------------------------------------------
SetSourceConnection(int id,vtkAlgorithmOutput * algOutput)495 void vtkTensorGlyph::SetSourceConnection(int id, vtkAlgorithmOutput* algOutput)
496 {
497 if (id < 0)
498 {
499 vtkErrorMacro("Bad index " << id << " for source.");
500 return;
501 }
502
503 int numConnections = this->GetNumberOfInputConnections(1);
504 if (id < numConnections)
505 {
506 this->SetNthInputConnection(1, id, algOutput);
507 }
508 else if (id == numConnections && algOutput)
509 {
510 this->AddInputConnection(1, algOutput);
511 }
512 else if (algOutput)
513 {
514 vtkWarningMacro("The source id provided is larger than the maximum "
515 "source id, using " << numConnections << " instead.");
516 this->AddInputConnection(1, algOutput);
517 }
518 }
519
520 //----------------------------------------------------------------------------
SetSourceData(vtkPolyData * source)521 void vtkTensorGlyph::SetSourceData(vtkPolyData *source)
522 {
523 this->SetInputData(1, source);
524 }
525
526 //----------------------------------------------------------------------------
GetSource()527 vtkPolyData *vtkTensorGlyph::GetSource()
528 {
529 if (this->GetNumberOfInputConnections(1) < 1)
530 {
531 return nullptr;
532 }
533 return vtkPolyData::SafeDownCast(this->GetExecutive()->GetInputData(1, 0));
534 }
535
536 //----------------------------------------------------------------------------
FillInputPortInformation(int port,vtkInformation * info)537 int vtkTensorGlyph::FillInputPortInformation(int port, vtkInformation *info)
538 {
539 if (port == 1)
540 {
541 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
542 return 1;
543 }
544 info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
545 return 1;
546 }
547
548 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)549 void vtkTensorGlyph::PrintSelf(ostream& os, vtkIndent indent)
550 {
551 this->Superclass::PrintSelf(os,indent);
552
553 os << indent << "Source: " << this->GetSource() << "\n";
554 os << indent << "Scaling: " << (this->Scaling ? "On\n" : "Off\n");
555 os << indent << "Scale Factor: " << this->ScaleFactor << "\n";
556 os << indent << "Extract Eigenvalues: " << (this->ExtractEigenvalues ? "On\n" : "Off\n");
557 os << indent << "Color Glyphs: " << (this->ColorGlyphs ? "On\n" : "Off\n");
558 os << indent << "Color Mode: " << this->ColorMode << endl;
559 os << indent << "Clamp Scaling: " << (this->ClampScaling ? "On\n" : "Off\n");
560 os << indent << "Max Scale Factor: " << this->MaxScaleFactor << "\n";
561 os << indent << "Three Glyphs: " << (this->ThreeGlyphs ? "On\n" : "Off\n");
562 os << indent << "Symmetric: " << (this->Symmetric ? "On\n" : "Off\n");
563 os << indent << "Length: " << this->Length << "\n";
564 }
565