You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
864 lines
23 KiB
864 lines
23 KiB
/*=========================================================================
|
|
|
|
Program: Visualization Toolkit
|
|
Module: $RCSfile: vtkHyperStreamline.cxx,v $
|
|
|
|
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
|
|
All rights reserved.
|
|
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
|
|
|
|
This software is distributed WITHOUT ANY WARRANTY; without even
|
|
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
|
PURPOSE. See the above copyright notice for more information.
|
|
|
|
=========================================================================*/
|
|
#include "vtkHyperStreamline.h"
|
|
|
|
#include "vtkCellArray.h"
|
|
#include "vtkDataSet.h"
|
|
#include "vtkFloatArray.h"
|
|
#include "vtkMath.h"
|
|
#include "vtkInformation.h"
|
|
#include "vtkInformationVector.h"
|
|
#include "vtkObjectFactory.h"
|
|
#include "vtkPointData.h"
|
|
#include "vtkPolyData.h"
|
|
|
|
vtkCxxRevisionMacro(vtkHyperStreamline, "$Revision: 1.60 $");
|
|
vtkStandardNewMacro(vtkHyperStreamline);
|
|
|
|
//
|
|
// Special classes for manipulating data
|
|
//BTX
|
|
class vtkHyperPoint { //;prevent man page generation
|
|
public:
|
|
vtkHyperPoint(); // method sets up storage
|
|
vtkHyperPoint &operator=(const vtkHyperPoint& hp); //for resizing
|
|
|
|
double X[3]; // position
|
|
vtkIdType CellId; // cell
|
|
int SubId; // cell sub id
|
|
double P[3]; // parametric coords in cell
|
|
double W[3]; // eigenvalues (sorted in decreasing value)
|
|
double *V[3]; // pointers to eigenvectors (also sorted)
|
|
double V0[3]; // storage for eigenvectors
|
|
double V1[3];
|
|
double V2[3];
|
|
double S; // scalar value
|
|
double D; // distance travelled so far
|
|
};
|
|
//ETX
|
|
|
|
class vtkHyperArray { //;prevent man page generation
|
|
public:
|
|
vtkHyperArray();
|
|
~vtkHyperArray()
|
|
{
|
|
if (this->Array)
|
|
{
|
|
delete [] this->Array;
|
|
}
|
|
};
|
|
vtkIdType GetNumberOfPoints() {return this->MaxId + 1;};
|
|
vtkHyperPoint *GetHyperPoint(vtkIdType i) {return this->Array + i;};
|
|
vtkHyperPoint *InsertNextHyperPoint()
|
|
{
|
|
if ( ++this->MaxId >= this->Size )
|
|
{
|
|
this->Resize(this->MaxId);
|
|
}
|
|
return this->Array + this->MaxId;
|
|
}
|
|
vtkHyperPoint *Resize(vtkIdType sz); //reallocates data
|
|
void Reset() {this->MaxId = -1;};
|
|
|
|
vtkHyperPoint *Array; // pointer to data
|
|
vtkIdType MaxId; // maximum index inserted thus far
|
|
vtkIdType Size; // allocated size of data
|
|
vtkIdType Extend; // grow array by this amount
|
|
double Direction; // integration direction
|
|
};
|
|
|
|
#define VTK_START_FROM_POSITION 0
|
|
#define VTK_START_FROM_LOCATION 1
|
|
|
|
vtkHyperPoint::vtkHyperPoint()
|
|
{
|
|
this->V[0] = this->V0;
|
|
this->V[1] = this->V1;
|
|
this->V[2] = this->V2;
|
|
}
|
|
|
|
vtkHyperPoint& vtkHyperPoint::operator=(const vtkHyperPoint& hp)
|
|
{
|
|
int i, j;
|
|
|
|
for (i=0; i<3; i++)
|
|
{
|
|
this->X[i] = hp.X[i];
|
|
this->P[i] = hp.P[i];
|
|
this->W[i] = hp.W[i];
|
|
for (j=0; j<3; j++)
|
|
{
|
|
this->V[j][i] = hp.V[j][i];
|
|
}
|
|
}
|
|
this->CellId = hp.CellId;
|
|
this->SubId = hp.SubId;
|
|
this->S = hp.S;
|
|
this->D = hp.D;
|
|
|
|
return *this;
|
|
}
|
|
|
|
vtkHyperArray::vtkHyperArray()
|
|
{
|
|
this->MaxId = -1;
|
|
this->Array = new vtkHyperPoint[1000];
|
|
this->Size = 1000;
|
|
this->Extend = 5000;
|
|
this->Direction = VTK_INTEGRATE_FORWARD;
|
|
}
|
|
|
|
vtkHyperPoint *vtkHyperArray::Resize(vtkIdType sz)
|
|
{
|
|
vtkHyperPoint *newArray;
|
|
vtkIdType newSize, i;
|
|
|
|
if (sz >= this->Size)
|
|
{
|
|
newSize = this->Size +
|
|
this->Extend*(((sz-this->Size)/this->Extend)+1);
|
|
}
|
|
else
|
|
{
|
|
newSize = sz;
|
|
}
|
|
|
|
newArray = new vtkHyperPoint[newSize];
|
|
|
|
for (i=0; i<sz; i++)
|
|
{
|
|
newArray[i] = this->Array[i];
|
|
}
|
|
|
|
this->Size = newSize;
|
|
delete [] this->Array;
|
|
this->Array = newArray;
|
|
|
|
return this->Array;
|
|
}
|
|
|
|
// Construct object with initial starting position (0,0,0); integration step
|
|
// length 0.2; step length 0.01; forward integration; terminal eigenvalue 0.0;
|
|
// number of sides 6; radius 0.5; and logarithmic scaling off.
|
|
vtkHyperStreamline::vtkHyperStreamline()
|
|
{
|
|
this->StartFrom = VTK_START_FROM_POSITION;
|
|
this->StartPosition[0] = this->StartPosition[1] = this->StartPosition[2] = 0.0;
|
|
|
|
this->StartCell = 0;
|
|
this->StartSubId = 0;
|
|
this->StartPCoords[0] = this->StartPCoords[1] = this->StartPCoords[2] = 0.5;
|
|
|
|
this->Streamers = NULL;
|
|
|
|
this->MaximumPropagationDistance = 100.0;
|
|
this->IntegrationStepLength = 0.2;
|
|
this->StepLength = 0.01;
|
|
this->IntegrationDirection = VTK_INTEGRATE_FORWARD;
|
|
this->TerminalEigenvalue = 0.0;
|
|
this->NumberOfSides = 6;
|
|
this->Radius = 0.5;
|
|
this->LogScaling = 0;
|
|
this->IntegrationEigenvector = VTK_INTEGRATE_MAJOR_EIGENVECTOR;
|
|
}
|
|
|
|
vtkHyperStreamline::~vtkHyperStreamline()
|
|
{
|
|
if ( this->Streamers )
|
|
{
|
|
delete [] this->Streamers;
|
|
}
|
|
}
|
|
|
|
// Specify the start of the hyperstreamline in the cell coordinate system.
|
|
// That is, cellId and subId (if composite cell), and parametric coordinates.
|
|
void vtkHyperStreamline::SetStartLocation(vtkIdType cellId, int subId,
|
|
double pcoords[3])
|
|
{
|
|
if ( cellId != this->StartCell || subId != this->StartSubId ||
|
|
pcoords[0] != this->StartPCoords[0] ||
|
|
pcoords[1] != this->StartPCoords[1] ||
|
|
pcoords[2] != this->StartPCoords[2] )
|
|
{
|
|
this->Modified();
|
|
this->StartFrom = VTK_START_FROM_LOCATION;
|
|
|
|
this->StartCell = cellId;
|
|
this->StartSubId = subId;
|
|
this->StartPCoords[0] = pcoords[0];
|
|
this->StartPCoords[1] = pcoords[1];
|
|
this->StartPCoords[2] = pcoords[2];
|
|
}
|
|
}
|
|
|
|
// Specify the start of the hyperstreamline in the cell coordinate system.
|
|
// That is, cellId and subId (if composite cell), and parametric coordinates.
|
|
void vtkHyperStreamline::SetStartLocation(vtkIdType cellId, int subId,
|
|
double r, double s, double t)
|
|
{
|
|
double pcoords[3];
|
|
pcoords[0] = r;
|
|
pcoords[1] = s;
|
|
pcoords[2] = t;
|
|
|
|
this->SetStartLocation(cellId, subId, pcoords);
|
|
}
|
|
|
|
// Get the starting location of the hyperstreamline in the cell coordinate
|
|
// system. Returns the cell that the starting point is in.
|
|
vtkIdType vtkHyperStreamline::GetStartLocation(int& subId, double pcoords[3])
|
|
{
|
|
subId = this->StartSubId;
|
|
pcoords[0] = this->StartPCoords[0];
|
|
pcoords[1] = this->StartPCoords[1];
|
|
pcoords[2] = this->StartPCoords[2];
|
|
return this->StartCell;
|
|
}
|
|
|
|
// Specify the start of the hyperstreamline in the global coordinate system.
|
|
// Starting from position implies that a search must be performed to find
|
|
// initial cell to start integration from.
|
|
void vtkHyperStreamline::SetStartPosition(double x[3])
|
|
{
|
|
if ( x[0] != this->StartPosition[0] || x[1] != this->StartPosition[1] ||
|
|
x[2] != this->StartPosition[2] )
|
|
{
|
|
this->Modified();
|
|
this->StartFrom = VTK_START_FROM_POSITION;
|
|
|
|
this->StartPosition[0] = x[0];
|
|
this->StartPosition[1] = x[1];
|
|
this->StartPosition[2] = x[2];
|
|
}
|
|
}
|
|
|
|
// Specify the start of the hyperstreamline in the global coordinate system.
|
|
// Starting from position implies that a search must be performed to find
|
|
// initial cell to start integration from.
|
|
void vtkHyperStreamline::SetStartPosition(double x, double y, double z)
|
|
{
|
|
double pos[3];
|
|
pos[0] = x;
|
|
pos[1] = y;
|
|
pos[2] = z;
|
|
|
|
this->SetStartPosition(pos);
|
|
}
|
|
|
|
// Get the start position of the hyperstreamline in global x-y-z coordinates.
|
|
double *vtkHyperStreamline::GetStartPosition()
|
|
{
|
|
return this->StartPosition;
|
|
}
|
|
|
|
// Make sure coordinate systems are consistent
|
|
static void FixVectors(double **prev, double **current, int iv, int ix, int iy)
|
|
{
|
|
double p0[3], p1[3], p2[3];
|
|
double v0[3], v1[3], v2[3];
|
|
double temp[3];
|
|
int i;
|
|
|
|
for (i=0; i<3; i++)
|
|
{
|
|
v0[i] = current[i][iv];
|
|
v1[i] = current[i][ix];
|
|
v2[i] = current[i][iy];
|
|
}
|
|
|
|
if ( prev == NULL ) //make sure coord system is right handed
|
|
{
|
|
vtkMath::Cross(v0,v1,temp);
|
|
if ( vtkMath::Dot(v2,temp) < 0.0 )
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
current[i][iy] *= -1.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
else //make sure vectors consistent from one point to the next
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
p0[i] = prev[i][iv];
|
|
p1[i] = prev[i][ix];
|
|
p2[i] = prev[i][iy];
|
|
}
|
|
if ( vtkMath::Dot(p0,v0) < 0.0 )
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
current[i][iv] *= -1.0;
|
|
}
|
|
}
|
|
if ( vtkMath::Dot(p1,v1) < 0.0 )
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
current[i][ix] *= -1.0;
|
|
}
|
|
}
|
|
if ( vtkMath::Dot(p2,v2) < 0.0 )
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
current[i][iy] *= -1.0;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
int vtkHyperStreamline::RequestData(
|
|
vtkInformation *vtkNotUsed(request),
|
|
vtkInformationVector **inputVector,
|
|
vtkInformationVector *outputVector)
|
|
{
|
|
// get the info objects
|
|
vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
|
|
vtkInformation *outInfo = outputVector->GetInformationObject(0);
|
|
|
|
// get the input and ouptut
|
|
vtkDataSet *input = vtkDataSet::SafeDownCast(
|
|
inInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
vtkPolyData *output = vtkPolyData::SafeDownCast(
|
|
outInfo->Get(vtkDataObject::DATA_OBJECT()));
|
|
|
|
vtkPointData *pd=input->GetPointData();
|
|
vtkDataArray *inScalars;
|
|
vtkDataArray *inTensors;
|
|
double tensor[9];
|
|
vtkHyperPoint *sNext, *sPtr;
|
|
int i, j, k, ptId, subId, iv, ix, iy;
|
|
vtkCell *cell;
|
|
double ev[3], xNext[3];
|
|
double d, step, dir, tol2, p[3];
|
|
double *w;
|
|
double dist2;
|
|
double closestPoint[3];
|
|
double *m[3], *v[3];
|
|
double m0[3], m1[3], m2[3];
|
|
double v0[3], v1[3], v2[3];
|
|
vtkDataArray *cellTensors;
|
|
vtkDataArray *cellScalars;
|
|
// set up working matrices
|
|
v[0] = v0; v[1] = v1; v[2] = v2;
|
|
m[0] = m0; m[1] = m1; m[2] = m2;
|
|
|
|
vtkDebugMacro(<<"Generating hyperstreamline(s)");
|
|
this->NumberOfStreamers = 0;
|
|
|
|
if ( ! (inTensors=pd->GetTensors()) )
|
|
{
|
|
vtkErrorMacro(<<"No tensor data defined!");
|
|
return 1;
|
|
}
|
|
w = new double[input->GetMaxCellSize()];
|
|
|
|
inScalars = pd->GetScalars();
|
|
|
|
cellTensors = vtkDataArray::CreateDataArray(inTensors->GetDataType());
|
|
cellScalars = vtkDataArray::CreateDataArray(inScalars->GetDataType());
|
|
int numComp;
|
|
if (inTensors)
|
|
{
|
|
numComp = inTensors->GetNumberOfComponents();
|
|
cellTensors->SetNumberOfComponents(numComp);
|
|
cellTensors->SetNumberOfTuples(VTK_CELL_SIZE);
|
|
}
|
|
if (inScalars)
|
|
{
|
|
numComp = inScalars->GetNumberOfComponents();
|
|
cellScalars->SetNumberOfComponents(numComp);
|
|
cellScalars->SetNumberOfTuples(VTK_CELL_SIZE);
|
|
}
|
|
|
|
tol2 = input->GetLength() / 1000.0;
|
|
tol2 = tol2 * tol2;
|
|
iv = this->IntegrationEigenvector;
|
|
ix = (iv + 1) % 3;
|
|
iy = (iv + 2) % 3;
|
|
//
|
|
// Create starting points
|
|
//
|
|
this->NumberOfStreamers = 1;
|
|
|
|
if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
|
|
{
|
|
this->NumberOfStreamers *= 2;
|
|
}
|
|
|
|
this->Streamers = new vtkHyperArray[this->NumberOfStreamers];
|
|
|
|
if ( this->StartFrom == VTK_START_FROM_POSITION )
|
|
{
|
|
sPtr = this->Streamers[0].InsertNextHyperPoint();
|
|
for (i=0; i<3; i++)
|
|
{
|
|
sPtr->X[i] = this->StartPosition[i];
|
|
}
|
|
sPtr->CellId = input->FindCell(this->StartPosition, NULL, (-1), 0.0,
|
|
sPtr->SubId, sPtr->P, w);
|
|
}
|
|
|
|
else //VTK_START_FROM_LOCATION
|
|
{
|
|
sPtr = this->Streamers[0].InsertNextHyperPoint();
|
|
cell = input->GetCell(sPtr->CellId);
|
|
cell->EvaluateLocation(sPtr->SubId, sPtr->P, sPtr->X, w);
|
|
}
|
|
//
|
|
// Finish initializing each hyperstreamline
|
|
//
|
|
this->Streamers[0].Direction = 1.0;
|
|
sPtr = this->Streamers[0].GetHyperPoint(0);
|
|
sPtr->D = 0.0;
|
|
if ( sPtr->CellId >= 0 ) //starting point in dataset
|
|
{
|
|
cell = input->GetCell(sPtr->CellId);
|
|
cell->EvaluateLocation(sPtr->SubId, sPtr->P, xNext, w);
|
|
|
|
inTensors->GetTuples(cell->PointIds, cellTensors);
|
|
|
|
// interpolate tensor, compute eigenfunctions
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] = 0.0;
|
|
}
|
|
}
|
|
for (k=0; k < cell->GetNumberOfPoints(); k++)
|
|
{
|
|
cellTensors->GetTuple(k, tensor);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] += tensor[i+3*j] * w[k];
|
|
}
|
|
}
|
|
}
|
|
|
|
vtkMath::Jacobi(m, sPtr->W, sPtr->V);
|
|
FixVectors(NULL, sPtr->V, iv, ix, iy);
|
|
|
|
if ( inScalars )
|
|
{
|
|
inScalars->GetTuples(cell->PointIds, cellScalars);
|
|
for (sPtr->S=0, i=0; i < cell->GetNumberOfPoints(); i++)
|
|
{
|
|
sPtr->S += cellScalars->GetTuple(i)[0] * w[i];
|
|
}
|
|
}
|
|
|
|
if ( this->IntegrationDirection == VTK_INTEGRATE_BOTH_DIRECTIONS )
|
|
{
|
|
this->Streamers[1].Direction = -1.0;
|
|
sNext = this->Streamers[1].InsertNextHyperPoint();
|
|
*sNext = *sPtr;
|
|
}
|
|
else if ( this->IntegrationDirection == VTK_INTEGRATE_BACKWARD )
|
|
{
|
|
this->Streamers[0].Direction = -1.0;
|
|
}
|
|
} //for hyperstreamline in dataset
|
|
//
|
|
// For each hyperstreamline, integrate in appropriate direction (using RK2).
|
|
//
|
|
for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
|
|
{
|
|
//get starting step
|
|
sPtr = this->Streamers[ptId].GetHyperPoint(0);
|
|
if ( sPtr->CellId < 0 )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
dir = this->Streamers[ptId].Direction;
|
|
cell = input->GetCell(sPtr->CellId);
|
|
cell->EvaluateLocation(sPtr->SubId, sPtr->P, xNext, w);
|
|
step = this->IntegrationStepLength * sqrt((double)cell->GetLength2());
|
|
inTensors->GetTuples(cell->PointIds, cellTensors);
|
|
if ( inScalars ) {inScalars->GetTuples(cell->PointIds, cellScalars);}
|
|
|
|
//integrate until distance has been exceeded
|
|
while ( sPtr->CellId >= 0 && fabs(sPtr->W[0]) > this->TerminalEigenvalue &&
|
|
sPtr->D < this->MaximumPropagationDistance )
|
|
{
|
|
|
|
//compute updated position using this step (Euler integration)
|
|
for (i=0; i<3; i++)
|
|
{
|
|
xNext[i] = sPtr->X[i] + dir * step * sPtr->V[i][iv];
|
|
}
|
|
|
|
//compute updated position using updated step
|
|
cell->EvaluatePosition(xNext, closestPoint, subId, p, dist2, w);
|
|
|
|
//interpolate tensor
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] = 0.0;
|
|
}
|
|
}
|
|
for (k=0; k < cell->GetNumberOfPoints(); k++)
|
|
{
|
|
cellTensors->GetTuple(k, tensor);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] += tensor[i+3*j] * w[k];
|
|
}
|
|
}
|
|
}
|
|
|
|
vtkMath::Jacobi(m, ev, v);
|
|
FixVectors(sPtr->V, v, iv, ix, iy);
|
|
|
|
//now compute final position
|
|
for (i=0; i<3; i++)
|
|
{
|
|
xNext[i] = sPtr->X[i] +
|
|
dir * (step/2.0) * (sPtr->V[i][iv] + v[i][iv]);
|
|
}
|
|
sNext = this->Streamers[ptId].InsertNextHyperPoint();
|
|
|
|
if ( cell->EvaluatePosition(xNext, closestPoint, sNext->SubId,
|
|
sNext->P, dist2, w) )
|
|
{ //integration still in cell
|
|
for (i=0; i<3; i++)
|
|
{
|
|
sNext->X[i] = closestPoint[i];
|
|
}
|
|
sNext->CellId = sPtr->CellId;
|
|
sNext->SubId = sPtr->SubId;
|
|
}
|
|
else
|
|
{ //integration has passed out of cell
|
|
sNext->CellId = input->FindCell(xNext, cell, sPtr->CellId, tol2,
|
|
sNext->SubId, sNext->P, w);
|
|
if ( sNext->CellId >= 0 ) //make sure not out of dataset
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
sNext->X[i] = xNext[i];
|
|
}
|
|
cell = input->GetCell(sNext->CellId);
|
|
inTensors->GetTuples(cell->PointIds, cellTensors);
|
|
if (inScalars){inScalars->GetTuples(cell->PointIds, cellScalars);}
|
|
step = this->IntegrationStepLength * sqrt((double)cell->GetLength2());
|
|
}
|
|
}
|
|
|
|
if ( sNext->CellId >= 0 )
|
|
{
|
|
cell->EvaluateLocation(sNext->SubId, sNext->P, xNext, w);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] = 0.0;
|
|
}
|
|
}
|
|
for (k=0; k < cell->GetNumberOfPoints(); k++)
|
|
{
|
|
cellTensors->GetTuple(k, tensor);
|
|
for (j=0; j<3; j++)
|
|
{
|
|
for (i=0; i<3; i++)
|
|
{
|
|
m[i][j] += tensor[i+3*j] * w[k];
|
|
}
|
|
}
|
|
}
|
|
|
|
vtkMath::Jacobi(m, sNext->W, sNext->V);
|
|
FixVectors(sPtr->V, sNext->V, iv, ix, iy);
|
|
|
|
if ( inScalars )
|
|
{
|
|
for (sNext->S=0.0, i=0; i < cell->GetNumberOfPoints(); i++)
|
|
{
|
|
sNext->S += cellScalars->GetTuple(i)[0] * w[i];
|
|
}
|
|
}
|
|
d = sqrt((double)vtkMath::Distance2BetweenPoints(sPtr->X,sNext->X));
|
|
sNext->D = sPtr->D + d;
|
|
}
|
|
|
|
sPtr = sNext;
|
|
|
|
}//for elapsed time
|
|
|
|
} //for each hyperstreamline
|
|
|
|
int retval = this->BuildTube(input, output);
|
|
|
|
delete [] w;
|
|
cellTensors->Delete();
|
|
cellScalars->Delete();
|
|
|
|
return retval;
|
|
}
|
|
|
|
int vtkHyperStreamline::BuildTube(vtkDataSet *input, vtkPolyData *output)
|
|
{
|
|
vtkHyperPoint *sPrev, *sPtr;
|
|
vtkPoints *newPts;
|
|
vtkFloatArray *newVectors;
|
|
vtkFloatArray *newNormals;
|
|
vtkFloatArray *newScalars=NULL;
|
|
vtkCellArray *newStrips;
|
|
vtkIdType i, npts, ptOffset=0;
|
|
int ptId, j, id, k, i1, i2;
|
|
double dOffset, x[3], v[3], s, r, r1[3], r2[3], stepLength;
|
|
double xT[3], sFactor, normal[3], w[3];
|
|
double theta=2.0*vtkMath::Pi()/this->NumberOfSides;
|
|
vtkPointData *outPD;
|
|
int iv, ix, iy;
|
|
vtkIdType numIntPts;
|
|
//
|
|
// Initialize
|
|
//
|
|
vtkDebugMacro(<<"Creating hyperstreamline tube");
|
|
if ( this->NumberOfStreamers <= 0 )
|
|
{
|
|
return 0;
|
|
}
|
|
|
|
stepLength = input->GetLength() * this->StepLength;
|
|
outPD = output->GetPointData();
|
|
|
|
iv = this->IntegrationEigenvector;
|
|
ix = (iv+1) % 3;
|
|
iy = (iv+2) % 3;
|
|
//
|
|
// Allocate
|
|
//
|
|
newPts = vtkPoints::New();
|
|
newPts ->Allocate(2500);
|
|
if ( input->GetPointData()->GetScalars() )
|
|
{
|
|
newScalars = vtkFloatArray::New();
|
|
newScalars->Allocate(2500);
|
|
}
|
|
newVectors = vtkFloatArray::New();
|
|
newVectors->SetNumberOfComponents(3);
|
|
newVectors->Allocate(7500);
|
|
newNormals = vtkFloatArray::New();
|
|
newNormals->SetNumberOfComponents(3);
|
|
newNormals->Allocate(7500);
|
|
newStrips = vtkCellArray::New();
|
|
newStrips->Allocate(newStrips->EstimateSize(3*this->NumberOfStreamers,
|
|
VTK_CELL_SIZE));
|
|
//
|
|
// Loop over all hyperstreamlines generating points
|
|
//
|
|
for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
|
|
{
|
|
if ( (numIntPts=this->Streamers[ptId].GetNumberOfPoints()) < 2 )
|
|
{
|
|
continue;
|
|
}
|
|
sPrev = this->Streamers[ptId].GetHyperPoint(0);
|
|
sPtr = this->Streamers[ptId].GetHyperPoint(1);
|
|
|
|
// compute scale factor
|
|
i = (sPrev->W[ix] > sPrev->W[iy] ? ix : iy);
|
|
if ( sPrev->W[i] == 0.0 )
|
|
{
|
|
sFactor = 1.0;
|
|
}
|
|
else
|
|
{
|
|
sFactor = this->Radius / sPrev->W[i];
|
|
}
|
|
|
|
if ( numIntPts == 2 && sPtr->CellId < 0 )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
dOffset = sPrev->D;
|
|
|
|
for ( npts=0, i=1; i < numIntPts && sPtr->CellId >= 0;
|
|
i++, sPrev=sPtr, sPtr=this->Streamers[ptId].GetHyperPoint(i) )
|
|
{
|
|
//
|
|
// Bracket steps and construct tube points
|
|
//
|
|
while ( dOffset >= sPrev->D && dOffset < sPtr->D )
|
|
{
|
|
r = (dOffset - sPrev->D) / (sPtr->D - sPrev->D);
|
|
|
|
for (j=0; j<3; j++) //compute point in center of tube
|
|
{
|
|
x[j] = sPrev->X[j] + r * (sPtr->X[j] - sPrev->X[j]);
|
|
v[j] = sPrev->V[j][iv] + r * (sPtr->V[j][iv] - sPrev->V[j][iv]);
|
|
r1[j] = sPrev->V[j][ix] + r * (sPtr->V[j][ix] - sPrev->V[j][ix]);
|
|
r2[j] = sPrev->V[j][iy] + r * (sPtr->V[j][iy] - sPrev->V[j][iy]);
|
|
w[j] = sPrev->W[j] + r * (sPtr->W[j] - sPrev->W[j]);
|
|
}
|
|
|
|
// construct points around tube
|
|
for (k=0; k < this->NumberOfSides; k++)
|
|
{
|
|
for (j=0; j<3; j++)
|
|
{
|
|
normal[j] = w[ix]*r1[j]*cos((double)k*theta) +
|
|
w[iy]*r2[j]*sin((double)k*theta);
|
|
xT[j] = x[j] + sFactor * normal[j];
|
|
}
|
|
id = newPts->InsertNextPoint(xT);
|
|
newVectors->InsertTuple(id,v);
|
|
vtkMath::Normalize(normal);
|
|
newNormals->InsertTuple(id,normal);
|
|
}
|
|
|
|
if ( newScalars ) //add scalars around tube
|
|
{
|
|
s = sPrev->S + r * (sPtr->S - sPrev->S);
|
|
for (k=0; k<this->NumberOfSides; k++)
|
|
{
|
|
newScalars->InsertNextTuple(&s);
|
|
}
|
|
}
|
|
|
|
npts++;
|
|
dOffset += stepLength;
|
|
|
|
} //while
|
|
} //for this hyperstreamline
|
|
|
|
//
|
|
// Generate the strips for this hyperstreamline
|
|
//
|
|
for (k=0; k<this->NumberOfSides; k++)
|
|
{
|
|
i1 = (k+1) % this->NumberOfSides;
|
|
newStrips->InsertNextCell(npts*2);
|
|
for (i=0; i < npts; i++)
|
|
{
|
|
//make sure strip definition consistent with normals
|
|
if (this->Streamers[ptId].Direction > 0.0)
|
|
{
|
|
i2 = i*this->NumberOfSides;
|
|
}
|
|
else
|
|
{
|
|
i2 = (npts - i - 1) * this->NumberOfSides;
|
|
}
|
|
newStrips->InsertCellPoint(ptOffset+i2+k);
|
|
newStrips->InsertCellPoint(ptOffset+i2+i1);
|
|
}
|
|
}//for all tube sides
|
|
|
|
ptOffset += this->NumberOfSides*npts;
|
|
|
|
} //for all hyperstreamlines
|
|
//
|
|
// Update ourselves
|
|
//
|
|
output->SetPoints(newPts);
|
|
newPts->Delete();
|
|
|
|
output->SetStrips(newStrips);
|
|
newStrips->Delete();
|
|
|
|
if ( newScalars )
|
|
{
|
|
int idx = outPD->AddArray(newScalars);
|
|
outPD->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
|
|
newScalars->Delete();
|
|
}
|
|
|
|
outPD->SetNormals(newNormals);
|
|
newNormals->Delete();
|
|
|
|
outPD->SetVectors(newVectors);
|
|
newVectors->Delete();
|
|
|
|
output->Squeeze();
|
|
|
|
return 1;
|
|
}
|
|
|
|
int vtkHyperStreamline::FillInputPortInformation(int, vtkInformation *info)
|
|
{
|
|
info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
|
|
return 1;
|
|
}
|
|
|
|
void vtkHyperStreamline::PrintSelf(ostream& os, vtkIndent indent)
|
|
{
|
|
this->Superclass::PrintSelf(os,indent);
|
|
|
|
if ( this->StartFrom == VTK_START_FROM_POSITION )
|
|
{
|
|
os << indent << "Starting Position: (" << this->StartPosition[0] << ","
|
|
<< this->StartPosition[1] << ", " << this->StartPosition[2] << ")\n";
|
|
}
|
|
else
|
|
{
|
|
os << indent << "Starting Location:\n\tCell: " << this->StartCell
|
|
<< "\n\tSubId: " << this->StartSubId << "\n\tP.Coordinates: ("
|
|
<< this->StartPCoords[0] << ", "
|
|
<< this->StartPCoords[1] << ", "
|
|
<< this->StartPCoords[2] << ")\n";
|
|
}
|
|
|
|
os << indent << "Maximum Propagation Distance: "
|
|
<< this->MaximumPropagationDistance << "\n";
|
|
|
|
if ( this->IntegrationDirection == VTK_INTEGRATE_FORWARD )
|
|
{
|
|
os << indent << "Integration Direction: FORWARD\n";
|
|
}
|
|
else if ( this->IntegrationDirection == VTK_INTEGRATE_BACKWARD )
|
|
{
|
|
os << indent << "Integration Direction: BACKWARD\n";
|
|
}
|
|
else
|
|
{
|
|
os << indent << "Integration Direction: FORWARD & BACKWARD\n";
|
|
}
|
|
|
|
os << indent << "Integration Step Length: " << this->IntegrationStepLength << "\n";
|
|
os << indent << "Step Length: " << this->StepLength << "\n";
|
|
|
|
os << indent << "Terminal Eigenvalue: " << this->TerminalEigenvalue << "\n";
|
|
|
|
os << indent << "Radius: " << this->Radius << "\n";
|
|
os << indent << "Number Of Sides: " << this->NumberOfSides << "\n";
|
|
os << indent << "Logarithmic Scaling: " << (this->LogScaling ? "On\n" : "Off\n");
|
|
|
|
if ( this->IntegrationEigenvector == 0 )
|
|
{
|
|
os << indent << "Integrate Along Major Eigenvector\n";
|
|
}
|
|
else if ( this->IntegrationEigenvector == 1 )
|
|
{
|
|
os << indent << "Integrate Along Medium Eigenvector\n";
|
|
}
|
|
else
|
|
{
|
|
os << indent << "Integrate Along Minor Eigenvector\n";
|
|
}
|
|
}
|
|
|