95{
96
97 vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
98 vtkInformation* sourceInfo = inputVector[1]->GetInformationObject(0);
99 vtkInformation* outInfo = outputVector->GetInformationObject(0);
100
101
102 vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
103 vtkPolyData* source = vtkPolyData::SafeDownCast(sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
104 vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
105
106 vtkDataArray* inTensors;
107 double tensor[9];
108 vtkDataArray* inScalars;
109 vtkIdType numPts, numSourcePts, numSourceCells, inPtId, i;
110 int j;
111 vtkPoints* sourcePts;
112 vtkDataArray* sourceNormals;
113 vtkCellArray *sourceCells, *cells;
114 vtkPoints* newPts;
115 vtkFloatArray* newScalars = nullptr;
116 vtkFloatArray* newNormals = nullptr;
117 double x[3], s;
118 vtkTransform* trans;
119 vtkCell* cell;
120 vtkIdList* cellPts;
121 int npts;
122 vtkIdType* pts;
123 vtkIdType ptIncr, cellId;
124 vtkIdType subIncr;
125 int numDirs, dir, eigen_dir, symmetric_dir;
126 vtkMatrix4x4* matrix;
127 double *m[3], w[3], *v[3];
128 double m0[3], m1[3], m2[3];
129 double v0[3], v1[3], v2[3];
130 double xv[3], yv[3], zv[3];
131 double maxScale;
132
134
135
136 m[0] = m0;
137 m[1] = m1;
138 m[2] = m2;
139 v[0] = v0;
140 v[1] = v1;
141 v[2] = v2;
142
143 vtkDebugMacro(<< "Generating tensor glyphs");
144
145 vtkPointData* outPD = output->GetPointData();
146 inTensors = this->GetInputArrayToProcess(0, inputVector);
147 inScalars = this->GetInputArrayToProcess(1, inputVector);
148 numPts = input->GetNumberOfPoints();
149
150 if ((inTensors == nullptr) || numPts < 1) {
151 vtkErrorMacro(<< "No data to glyph!");
152 return 1;
153 }
154
155 pts = new vtkIdType[source->GetMaxCellSize()];
156 trans = vtkTransform::New();
157 matrix = vtkMatrix4x4::New();
158
159
160
161
162 sourcePts = source->GetPoints();
163 numSourcePts = sourcePts->GetNumberOfPoints();
164 numSourceCells = source->GetNumberOfCells();
165
166 newPts = vtkPoints::New();
167 newPts->Allocate(numDirs * numPts * numSourcePts);
168
169
170 if ((sourceCells = source->GetVerts())->GetNumberOfCells() > 0) {
171 cells = vtkCellArray::New();
172 cells->Allocate(numDirs * numPts * sourceCells->GetSize());
173 output->SetVerts(cells);
174 cells->Delete();
175 }
176 if ((sourceCells = this->
GetSource()->GetLines())->GetNumberOfCells() > 0) {
177 cells = vtkCellArray::New();
178 cells->Allocate(numDirs * numPts * sourceCells->GetSize());
179 output->SetLines(cells);
180 cells->Delete();
181 }
182 if ((sourceCells = this->
GetSource()->GetPolys())->GetNumberOfCells() > 0) {
183 cells = vtkCellArray::New();
184 cells->Allocate(numDirs * numPts * sourceCells->GetSize());
185 output->SetPolys(cells);
186 cells->Delete();
187 }
188 if ((sourceCells = this->
GetSource()->GetStrips())->GetNumberOfCells() > 0) {
189 cells = vtkCellArray::New();
190 cells->Allocate(numDirs * numPts * sourceCells->GetSize());
191 output->SetStrips(cells);
192 cells->Delete();
193 }
194
195
196 vtkPointData* pd = this->
GetSource()->GetPointData();
197
201 {
202 newScalars = vtkFloatArray::New();
203 newScalars->SetNumberOfComponents(4);
204 newScalars->Allocate(numDirs * numPts * numSourcePts);
206 newScalars->SetName("MaxEigenvalue");
207 }
208 else {
209 newScalars->SetName(inScalars->GetName());
210 }
211 }
212 else {
213 outPD->CopyAllOff();
214 outPD->CopyScalarsOn();
215 outPD->CopyAllocate(pd, numDirs * numPts * numSourcePts);
216 }
217 if ((sourceNormals = pd->GetNormals()) != nullptr) {
218 newNormals = vtkFloatArray::New();
219 newNormals->SetNumberOfComponents(3);
220 newNormals->SetName("Normals");
221 newNormals->Allocate(numDirs * 3 * numPts * numSourcePts);
222 }
223
224
225
226 for (inPtId = 0; inPtId < numPts; inPtId++) {
227 ptIncr = numDirs * inPtId * numSourcePts;
228 for (cellId = 0; cellId < numSourceCells; cellId++) {
229 cell = this->
GetSource()->GetCell(cellId);
230 cellPts = cell->GetPointIds();
231 npts = (int)cellPts->GetNumberOfIds();
232 for (dir = 0; dir < numDirs; dir++) {
233
234
235 subIncr = ptIncr + dir * numSourcePts;
236 for (i = 0; i < npts; i++) {
237 pts[i] = cellPts->GetId(i) + subIncr;
238 }
239 output->InsertNextCell(cell->GetCellType(), npts, pts);
240 }
241 }
242 }
243
244
245
246 trans->PreMultiply();
247
248 for (inPtId = 0; inPtId < numPts; inPtId++) {
249 ptIncr = numDirs * inPtId * numSourcePts;
250
251
252
253 inTensors->GetTuple(inPtId, tensor);
254 if (inTensors->GetNumberOfComponents() == 6) {
255 vtkMath::TensorFromSymmetricTensor(tensor);
256 }
257
258
260 {
261
262
263 for (j = 0; j < 3; j++) {
264 for (i = 0; i < 3; i++) {
265 m[i][j] = 0.5 * (tensor[i + 3 * j] + tensor[j + 3 * i]);
266 }
267 }
268 vtkMath::Jacobi(m, w, v);
269
270
271 xv[0] = v[0][0];
272 xv[1] = v[1][0];
273 xv[2] = v[2][0];
274 yv[0] = v[0][1];
275 yv[1] = v[1][1];
276 yv[2] = v[2][1];
277 zv[0] = v[0][2];
278 zv[1] = v[1][2];
279 zv[2] = v[2][2];
280 }
281 else
282 {
283 for (i = 0; i < 3; i++) {
284 xv[i] = tensor[i];
285 yv[i] = tensor[i + 3];
286 zv[i] = tensor[i + 6];
287 }
288 w[0] = vtkMath::Normalize(xv);
289 w[1] = vtkMath::Normalize(yv);
290 w[2] = vtkMath::Normalize(zv);
291 }
292
293
297
299 for (maxScale = 0.0, i = 0; i < 3; i++) {
300 if (maxScale < fabs(w[i])) {
301 maxScale = fabs(w[i]);
302 }
303 }
306 for (i = 0; i < 3; i++) {
307 w[i] *= maxScale;
308 }
309 }
310 }
311
312
313
314
315 for (maxScale = 0.0, i = 0; i < 3; i++) {
316 if (w[i] > maxScale) {
317 maxScale = w[i];
318 }
319 }
320 if (maxScale == 0.0) {
321 maxScale = 1.0;
322 }
323 for (i = 0; i < 3; i++) {
324 if (w[i] == 0.0) {
325 w[i] = maxScale * 1.0e-06;
326 }
327 }
328
329
330
331 for (dir = 0; dir < numDirs; dir++) {
332 eigen_dir = dir % (this->
ThreeGlyphs != 0 ? 3 : 1);
333 symmetric_dir = dir / (this->
ThreeGlyphs != 0 ? 3 : 1);
334
335
336 trans->Identity();
337
338
339 input->GetPoint(inPtId, x);
340 trans->Translate(x[0], x[1], x[2]);
341
342
343 matrix->Element[0][0] = xv[0];
344 matrix->Element[0][1] = yv[0];
345 matrix->Element[0][2] = zv[0];
346 matrix->Element[1][0] = xv[1];
347 matrix->Element[1][1] = yv[1];
348 matrix->Element[1][2] = zv[1];
349 matrix->Element[2][0] = xv[2];
350 matrix->Element[2][1] = yv[2];
351 matrix->Element[2][2] = zv[2];
352 trans->Concatenate(matrix);
353
354 if (eigen_dir == 1) {
355 trans->RotateZ(90.0);
356 }
357
358 if (eigen_dir == 2) {
359 trans->RotateY(-90.0);
360 }
361
364 }
365 else {
366 trans->Scale(w[0], w[1], w[2]);
367 }
368
369
370 if (symmetric_dir == 1) {
371 trans->Scale(-1., 1., 1.);
372 }
373
374
375
376
377
378 if (w[eigen_dir] < 0 && numDirs > 1) {
379 trans->Translate(-this->
Length, 0., 0.);
380 }
381
382
383
384 trans->TransformPoints(sourcePts, newPts);
385
386
387
388 if (newNormals != nullptr) {
389
390
391
392
393 if (trans->GetMatrix()->Determinant() < 0) {
394 trans->Scale(-1.0, -1.0, -1.0);
395 }
396 trans->TransformNormals(sourceNormals, newNormals);
397 }
398
399
400 if ((this->
ColorGlyphs != 0) && (inScalars !=
nullptr)
402 {
403 for (i = 0; i < numSourcePts; i++) {
404 auto st = inScalars->GetTuple(inPtId);
405 newScalars->InsertTuple4(ptIncr + i, st[0], st[1], st[2], st[3]);
406 }
407 }
409
410
411 s = w[eigen_dir];
412 for (i = 0; i < numSourcePts; i++) {
413 newScalars->InsertTuple(ptIncr + i, &s);
414 }
415 }
416 else {
417 for (i = 0; i < numSourcePts; i++) {
418 outPD->CopyData(pd, i, ptIncr + i);
419 }
420 }
421 ptIncr += numSourcePts;
422 }
423 }
424 vtkDebugMacro(<< "Generated " << numPts << " tensor glyphs");
425
426
427
428 delete[] pts;
429
430 output->SetPoints(newPts);
431 newPts->Delete();
432
433 if (newScalars != nullptr) {
434 int idx = outPD->AddArray(newScalars);
435 outPD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
436 newScalars->Delete();
437 }
438
439 if (newNormals != nullptr) {
440 outPD->SetNormals(newNormals);
441 newNormals->Delete();
442 }
443
444 output->Squeeze();
445 trans->Delete();
446 matrix->Delete();
447
448 return 1;
449}
vtkPolyData * GetSource()
vtkTypeBool ExtractEigenvalues