101{
102
103 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
104 vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
105 vtkInformation *outInfo = outputVector->GetInformationObject(0);
106
107
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
143
144
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
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
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
206 vtkPointData *pd = this->
GetSource()->GetPointData();
207
211 {
212 newScalars = vtkFloatArray::New();
213 newScalars->SetNumberOfComponents(4);
214 newScalars->Allocate(numDirs*numPts*numSourcePts);
216 {
217 newScalars->SetName("MaxEigenvalue");
218 }
219 else
220 {
221 newScalars->SetName(inScalars->GetName());
222 }
223 }
224 else
225 {
226 outPD->CopyAllOff();
227 outPD->CopyScalarsOn();
228 outPD->CopyAllocate(pd,numDirs*numPts*numSourcePts);
229 }
230 if ( (sourceNormals = pd->GetNormals()) )
231 {
232 newNormals = vtkFloatArray::New();
233 newNormals->SetNumberOfComponents(3);
234 newNormals->SetName("Normals");
235 newNormals->Allocate(numDirs*3*numPts*numSourcePts);
236 }
237
238
239
240 for (inPtId=0; inPtId < numPts; inPtId++)
241 {
242 ptIncr = numDirs * inPtId * numSourcePts;
243 for (cellId=0; cellId < numSourceCells; cellId++)
244 {
245 cell = this->
GetSource()->GetCell(cellId);
246 cellPts = cell->GetPointIds();
247 npts = (int)cellPts->GetNumberOfIds();
248 for (dir=0; dir < numDirs; dir++)
249 {
250
251
252 subIncr = ptIncr + dir*numSourcePts;
253 for (i=0; i < npts; i++)
254 {
255 pts[i] = cellPts->GetId(i) + subIncr;
256 }
257 output->InsertNextCell(cell->GetCellType(),npts,pts);
258 }
259 }
260 }
261
262
263
264 trans->PreMultiply();
265
266 for (inPtId=0; inPtId < numPts; inPtId++)
267 {
268 ptIncr = numDirs * inPtId * numSourcePts;
269
270
271
272 inTensors->GetTuple(inPtId, tensor);
273 if (inTensors->GetNumberOfComponents() == 6)
274 {
275 vtkMath::TensorFromSymmetricTensor(tensor);
276 }
277
278
280 {
281
282
283 for (j=0; j<3; j++)
284 {
285 for (i=0; i<3; i++)
286 {
287 m[i][j] = 0.5 * (tensor[i + 3 * j] + tensor[j + 3 * i]);
288 }
289 }
290 vtkMath::Jacobi(m, w, v);
291
292
293 xv[0] = v[0][0]; xv[1] = v[1][0]; xv[2] = v[2][0];
294 yv[0] = v[0][1]; yv[1] = v[1][1]; yv[2] = v[2][1];
295 zv[0] = v[0][2]; zv[1] = v[1][2]; zv[2] = v[2][2];
296 }
297 else
298 {
299 for (i=0; i<3; i++)
300 {
301 xv[i] = tensor[i];
302 yv[i] = tensor[i+3];
303 zv[i] = tensor[i+6];
304 }
305 w[0] = vtkMath::Normalize(xv);
306 w[1] = vtkMath::Normalize(yv);
307 w[2] = vtkMath::Normalize(zv);
308 }
309
310
314
316 {
317 for (maxScale=0.0, i=0; i<3; i++)
318 {
319 if ( maxScale < fabs(w[i]) )
320 {
321 maxScale = fabs(w[i]);
322 }
323 }
325 {
327 for (i=0; i<3; i++)
328 {
329 w[i] *= maxScale;
330 }
331 }
332 }
333
334
335
336
337 for (maxScale=0.0, i=0; i<3; i++)
338 {
339 if ( w[i] > maxScale )
340 {
341 maxScale = w[i];
342 }
343 }
344 if ( maxScale == 0.0 )
345 {
346 maxScale = 1.0;
347 }
348 for (i=0; i<3; i++)
349 {
350 if ( w[i] == 0.0 )
351 {
352 w[i] = maxScale * 1.0e-06;
353 }
354 }
355
356
357
358 for (dir=0; dir < numDirs; dir++)
359 {
362
363
364 trans->Identity();
365
366
367 input->GetPoint(inPtId, x);
368 trans->Translate(x[0], x[1], x[2]);
369
370
371 matrix->Element[0][0] = xv[0];
372 matrix->Element[0][1] = yv[0];
373 matrix->Element[0][2] = zv[0];
374 matrix->Element[1][0] = xv[1];
375 matrix->Element[1][1] = yv[1];
376 matrix->Element[1][2] = zv[1];
377 matrix->Element[2][0] = xv[2];
378 matrix->Element[2][1] = yv[2];
379 matrix->Element[2][2] = zv[2];
380 trans->Concatenate(matrix);
381
382 if (eigen_dir == 1)
383 {
384 trans->RotateZ(90.0);
385 }
386
387 if (eigen_dir == 2)
388 {
389 trans->RotateY(-90.0);
390 }
391
393 {
395 }
396 else
397 {
398 trans->Scale(w[0], w[1], w[2]);
399 }
400
401
402 if (symmetric_dir == 1)
403 {
404 trans->Scale(-1.,1.,1.);
405 }
406
407
408
409
410
411 if (w[eigen_dir] < 0 && numDirs > 1)
412 {
413 trans->Translate(-this->
Length, 0., 0.);
414 }
415
416
417
418 trans->TransformPoints(sourcePts,newPts);
419
420
421
422 if ( newNormals )
423 {
424
425
426
427
428 if (trans->GetMatrix()->Determinant() < 0)
429 {
430 trans->Scale(-1.0,-1.0,-1.0);
431 }
432 trans->TransformNormals(sourceNormals,newNormals);
433 }
434
435
438 {
439 for (i=0; i < numSourcePts; i++)
440 {
441 auto st = inScalars->GetTuple(inPtId);
442 newScalars->InsertTuple4(ptIncr+i,st[0],st[1],st[2],st[3]);
443 }
444 }
447 {
448
449
450 s = w[eigen_dir];
451 for (i=0; i < numSourcePts; i++)
452 {
453 newScalars->InsertTuple(ptIncr+i, &s);
454 }
455 }
456 else
457 {
458 for (i=0; i < numSourcePts; i++)
459 {
460 outPD->CopyData(pd,i,ptIncr+i);
461 }
462 }
463 ptIncr += numSourcePts;
464 }
465 }
466 vtkDebugMacro(<<"Generated " << numPts <<" tensor glyphs");
467
468
469
470 delete [] pts;
471
472 output->SetPoints(newPts);
473 newPts->Delete();
474
475 if ( newScalars )
476 {
477 int idx = outPD->AddArray(newScalars);
478 outPD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
479 newScalars->Delete();
480 }
481
482 if ( newNormals )
483 {
484 outPD->SetNormals(newNormals);
485 newNormals->Delete();
486 }
487
488 output->Squeeze();
489 trans->Delete();
490 matrix->Delete();
491
492 return 1;
493}
vtkPolyData * GetSource()
vtkTypeBool ExtractEigenvalues