[vtkusers] End locations of vtkStreamTracer
Tijmen Klein
T.R.Klein at student.rug.nl
Mon Oct 3 09:21:35 EDT 2011
Hello Rocco,
Thank you very much for your example. My question was still relevant, and I
was able to modify your example to my needs. I'll post my code here, just in
case someone else needs it too (perhaps in the future). It is based on
Rocco's code, and calculates the endpoints of a vtkStreamTracer.
vtkSmartPointer<vtkPoints>
getStreamTracerEndPoints(vtkSmartPointer<vtkStreamTracer> streamTracer)
{
vtkPolyData* inputStreamlineData = streamTracer->GetOutput();
inputStreamlineData->BuildCells();
inputStreamlineData->BuildLinks();
vtkSmartPointer<vtkPoints> endPoints =
vtkSmartPointer<vtkPoints>::New();
vtkSmartPointer<vtkIdList> currentLinePointIds =
vtkSmartPointer<vtkIdList>::New();
vtkSmartPointer<vtkIdList> incidentStreamlineIdsList =
vtkSmartPointer<vtkIdList>::New();
vtkCellArray *unprocessedStreamlines = inputStreamlineData->GetLines();
if (unprocessedStreamlines == NULL)
{
cerr << "No streamlines available ==> Can not get end points!" << endl;
} else
{
unprocessedStreamlines->InitTraversal();
vtkIdType npts, *pts;
while(unprocessedStreamlines->GetNextCell(npts, pts) && npts > 0)
{
incidentStreamlineIdsList->Reset();
currentLinePointIds->Reset();
inputStreamlineData->GetPointCells(pts[0], incidentStreamlineIdsList);
inputStreamlineData->GetCellPoints(incidentStreamlineIdsList->GetId(0),
currentLinePointIds);
int lastPointId =
currentLinePointIds->GetId(currentLinePointIds->GetNumberOfIds()-1);
endPoints->InsertNextPoint(inputStreamlineData->GetPoint(lastPointId));
}
}
return endPoints;
}
A small example on how to use it. Get the end points of a vtkStreamTracer
and use those points as the need seed points:
vtkSmartPointer<vtkPolyData> sourceData =
vtkSmartPointer<vtkPolyData>::New();
sourceData->SetPoints(getStreamTracerEndPoints(streamTracer));
streamTracer->SetSource(sourceData);
Cheers,
Tijmen
On Fri, Sep 30, 2011 at 9:39 AM, Rocco Gasteiger <post at rocco-gasteiger.de>wrote:
> Hello Tijmen,****
>
> ** **
>
> I do not know if you question is still relevant but I will give one
> suggestion how you can find the end points of your streamlines. The output
> of vtkStreamTracer is a vtkPolyDataSet which contains the streamlines as
> multiple line cells. You can iterate over each line to retrieve the vertices
> of that line. These vertices are ordered beginning with the first vertex in
> that line. ****
>
> ** **
>
> However, you must be careful! If you propagate in both directions within
> the vtkStreamTracer you have up to two lines for a "whole" streamline. One
> for the forward and one for the backward propagation. These two line
> segments are not connected! I wrote a function which starts with the first
> vertex of the current line. I checked this vertex for two cases. First case:
> this vertex is the start point of a streamline with only one segment, that
> means only forward or backward propagation was done. Second case: like the
> first case but with another vertex at the same position but it is the start
> point of the second line segment. I checked this by searching in the
> neighborhood of that vertex, with radius 0.0, if another vertex is there.
> For this searching I used vtkPointLocator. I connect the two line segments
> later in the right order and discard one of the two vertices. Below I give
> you the C++ code for that function. I hope that the comments are sufficient.
> If not do not hesitate to ask me.****
>
> ** **
>
> Best regards, Rocco****
>
> ** **
>
> void UMDStreamlineStyleMapper::sortFlowVertices()****
>
> {****
>
> inputStreamlineData->BuildCells();****
>
> inputStreamlineData->BuildLinks();****
>
> ** **
>
> vtkSmartPointer<vtkPolyData> sortedStreamlineData =
> vtkSmartPointer<vtkPolyData>::New();****
>
> vtkSmartPointer<vtkPoints> sortedStreamlinePoints =
> vtkSmartPointer<vtkPoints>::New();****
>
> vtkSmartPointer<vtkCellArray> sortedStreamlines =
> vtkSmartPointer<vtkCellArray>::New();****
>
> vtkSmartPointer<vtkDoubleArray> streamlineScalars =
> vtkSmartPointer<vtkDoubleArray>::New();****
>
> vtkSmartPointer<vtkPointLocator> pointLocator = vtkSmartPointer<
> vtkPointLocator>::New();****
>
> pointLocator->SetDataSet(inputStreamlineData);****
>
> pointLocator->BuildLocator();****
>
> vtkSmartPointer<vtkIdList> composedLinePointIDs =
> vtkSmartPointer<vtkIdList>::New();****
>
> vtkSmartPointer<vtkIdList> rearrangedLinePointIDs =
> vtkSmartPointer<vtkIdList>::New(); ****
>
> vtkSmartPointer<vtkIdList> identicalSeedPointIdsList =
> vtkSmartPointer<vtkIdList>::New();****
>
> vtkSmartPointer<vtkIdList> currentLinePointIds = vtkSmartPointer
> <vtkIdList>::New();****
>
> vtkSmartPointer<vtkIdList> incidentStreamlineIdsList =
> vtkSmartPointer<vtkIdList>::New();****
>
> vtkSmartPointer<vtkIdList> processedSeedpointIds =
> vtkSmartPointer<vtkIdList>::New();****
>
> ** **
>
> progressBar->setFloatValue(0.0);****
>
> ** **
>
> bool hasScalars = false;****
>
> float maxLength = 0.0;****
>
> int newPointId = -1;****
>
> int numProcessedStreamLines = -1;****
>
> ** **
>
> if (inputStreamlineData->GetPointData()->GetScalars())****
>
> {****
>
> streamlineScalars->SetName(inputStreamlineData->GetPointData
> ()->GetScalars()->GetName());****
>
> hasScalars = true;****
>
> } ****
>
> ** **
>
> vtkCellArray *unprocessedStreamlines = inputStreamlineData->GetLines
> ();****
>
> if (unprocessedStreamlines == NULL)****
>
> {****
>
> cerr << "No streamlines available ==> no sorting is
> performed!" << endl;****
>
> } else****
>
> {****
>
> unprocessedStreamlines->InitTraversal();****
>
> vtkIdType npts, *pts;****
>
> while(unprocessedStreamlines->GetNextCell(npts, pts) && npts >
> 0)****
>
> {****
>
> identicalSeedPointIdsList->Reset();****
>
> incidentStreamlineIdsList->Reset();****
>
> currentLinePointIds->Reset();****
>
> rearrangedLinePointIDs->Reset();****
>
> composedLinePointIDs->Reset();****
>
> ** **
>
> numProcessedStreamLines++;****
>
> // The first point of a streamline corresponds to its
> seed point. We store its id to skip processing a second streamline ****
>
> // propagation because we already have processed this
> line (see below)****
>
> int id = pts[0];****
>
> if (processedSeedpointIds->IsId(id) != -1)****
>
> {****
>
> // Continue with the next streamline starting
> from a new seed point****
>
> continue;****
>
> }****
>
> // An unprocessed seed point is found. Test if this
> point exists twice at the same position. ****
>
> // The second point is used as seed point for another
> streamline (forward or backward propagation). ****
>
> processedSeedpointIds->InsertNextId(id);****
>
> // Get the position of the current seed point****
>
> double seedPoint[3];****
>
> inputStreamlineData->GetPoint(id, seedPoint);****
>
> pointLocator->FindPointsWithinRadius(0.0, seedPoint,
> identicalSeedPointIdsList);****
>
> int numIdenticalSeedPoints = identicalSeedPointIdsList->
> GetNumberOfIds();****
>
> ****
>
> int secondId = 0;****
>
> int firstId = 0;****
>
> switch (numIdenticalSeedPoints)****
>
> {****
>
> case 1:****
>
> // Get all the point ids of the single streamline
> (in "identicalSeedPointIdsList" and "incidentStreamlineIdsList"****
>
> // are only one entry).****
>
> firstId = identicalSeedPointIdsList->GetId(0);****
>
> inputStreamlineData->GetPointCells(firstId,
> incidentStreamlineIdsList);****
>
> inputStreamlineData->GetCellPoints(
> incidentStreamlineIdsList->GetId(0), currentLinePointIds);****
>
> // Reverse streamline vertex ids because they
> start from the seed point.****
>
> for (int pointId = 0; pointId <
> currentLinePointIds->GetNumberOfIds(); pointId++)****
>
> {****
>
> // Reverse the ids because they start at the
> seed point.****
>
> int currentPointId = currentLinePointIds->
> GetId(currentLinePointIds->GetNumberOfIds()-pointId-1);****
>
> rearrangedLinePointIDs->InsertNextId(
> currentPointId);****
>
> }****
>
> break;****
>
> case 2:****
>
> // Get all the point ids of the first streamline
> (in "identicalSeedPointIdsList" are two entries****
>
> // and in "incidentStreamlineIdsList" is still one
> entry).****
>
> inputStreamlineData->GetPointCells(
> identicalSeedPointIdsList->GetId(0), incidentStreamlineIdsList);****
>
> inputStreamlineData->GetCellPoints(
> incidentStreamlineIdsList->GetId(0), currentLinePointIds);****
>
> // Reverse streamline vertex ids because they
> start from the seed point.****
>
> for (int pointId = 0; pointId <
> currentLinePointIds->GetNumberOfIds(); pointId++)****
>
> {****
>
> // Reverse the ids because they start at the
> seed point.****
>
> int currentPointId = currentLinePointIds->
> GetId(currentLinePointIds->GetNumberOfIds()-pointId-1);****
>
> rearrangedLinePointIDs->InsertNextId(
> currentPointId);****
>
> }****
>
> // Get all the point ids of the second streamline*
> ***
>
> currentLinePointIds->Reset();****
>
> incidentStreamlineIdsList->Reset();****
>
> inputStreamlineData->GetPointCells(
> identicalSeedPointIdsList->GetId(1), incidentStreamlineIdsList);****
>
> secondId = incidentStreamlineIdsList->GetId(0);***
> *
>
> inputStreamlineData->GetCellPoints(
> incidentStreamlineIdsList->GetId(0), currentLinePointIds);****
>
> ****
>
> // We do not reverse the vertex order because we
> continue the streamline direction further.****
>
> // We also skip the first point of the second
> streamline because this point is already inserted****
>
> // in the rearranged list.****
>
> for (int pointId = 1; pointId <
> currentLinePointIds->GetNumberOfIds(); pointId++)****
>
> {****
>
> // Reverse the ids because they start at the
> seed point.****
>
> int currentPointId = currentLinePointIds->
> GetId(pointId);****
>
> rearrangedLinePointIDs->InsertNextId(
> currentPointId);****
>
> } ****
>
> break;****
>
> default:****
>
> cerr << "UMDStreamlineStyleMapper::sortFlowVertices()==>
> Invalid number of streamlines at current seed point!" << endl;****
>
> break;****
>
> }****
>
> ** **
>
> float currentLength = 0.0; ****
>
> ** **
>
> for (int pointId = 0; pointId < rearrangedLinePointIDs->
> GetNumberOfIds(); pointId++)****
>
> {****
>
> newPointId++; ****
>
> ** **
>
> int currentPointId = rearrangedLinePointIDs->GetId
> (pointId);****
>
> ** **
>
> if (hasScalars)****
>
> {****
>
> double scalarValue = inputStreamlineData->
> GetPointData()->GetScalars()->GetTuple1(currentPointId);****
>
> streamlineScalars->InsertNextTuple1(
> scalarValue); ****
>
> }****
>
> ** **
>
> composedLinePointIDs->InsertNextId(newPointId);***
> *
>
> sortedStreamlinePoints->InsertNextPoint(
> inputStreamlineData->GetPoint(currentPointId));****
>
> ** **
>
> if (newPointId > 0)****
>
> {****
>
> double p0[3];****
>
> double p1[3]; ****
>
> sortedStreamlinePoints->GetPoint(newPointId-1,
> p0);****
>
> sortedStreamlinePoints->GetPoint(newPointId,
> p1);****
>
> currentLength += sqrt(vtkMath::
> Distance2BetweenPoints(p0, p1));****
>
> } ****
>
> } ****
>
> sortedStreamlines->InsertNextCell(composedLinePointIDs->
> GetNumberOfIds(), composedLinePointIDs->GetPointer(0));****
>
> ** **
>
> ** **
>
> if (currentLength > maxLength)****
>
> {****
>
> maxLength = currentLength;****
>
> }****
>
> setMaximumStreamlineLength(maxLength);****
>
> ** **
>
> // The first half of the progress state consists of
> sorting the flow vertices.****
>
> // The second half will be done during the generating of
> the edge mesh data structure.****
>
> progressBar->setFloatValue((float)
> numProcessedStreamLines/inputStreamlineData->GetNumberOfLines()); **
> **
>
> } // end while-loop****
>
> ** **
>
> cout << "sortedPoints->GetNumberOfPoints(): " <<
> sortedStreamlinePoints->GetNumberOfPoints() << endl;****
>
> //cout << "unsortedPolyLines->GetNumberOfPoints(): " <<
> unsortedPolyLines->GetNumberOfPoints() << endl;****
>
> //cout << "polyLines->GetNumberOfCells(): " <<
> polyLines->GetNumberOfCells() << endl;****
>
> //cout << "maxLength: " << maxLength << endl;****
>
> ** **
>
> sortedStreamlineData->SetPoints(sortedStreamlinePoints);****
>
> sortedStreamlineData->SetLines(sortedStreamlines);****
>
> sortedStreamlineData->GetPointData()->SetScalars(
> streamlineScalars);****
>
> sortedStreamlineData->Update();****
>
> inputStreamlineData = sortedStreamlineData;****
>
> vtkOutputSortedStreamlines.setNewOutputBaseFieldObject(
> sortedStreamlineData);****
>
> } ****
>
> }****
>
> ** **
>
> *Von:* vtkusers-bounces at vtk.org [mailto:vtkusers-bounces at vtk.org] *Im
> Auftrag von *Tijmen Klein
> *Gesendet:* Dienstag, 27. September 2011 15:10
> *An:* vtkusers at vtk.org
> *Betreff:* Re: [vtkusers] End locations of vtkStreamTracer****
>
> ** **
>
> A little update on my own question. I am now able to find the end location
> of a streamline, as long as I use 1 input point (setStartPosition() instead
> of setSource() ). I can find the end point using this trick:****
>
> ** **
>
> double* stopPosition;****
>
> long numPoints = streamer->GetOutput()->GetNumberOfPoints();****
>
> stopPosition = streamer->GetOutput()->GetPoint(numPoints-1);****
>
> ** **
>
> I can then use this stopPosition for the setStartPosition() method.
> Unfortunately, this does not work with multiple inputs. Is there anyone who
> could help me on that? Ofcouse, I could create an array of vtkStreamTracers,
> each one only tracing a single input point. This should work, but feels
> really ugly.****
>
> ** **
>
> Cheers,****
>
> Tijmen****
>
> ** **
>
> On Wed, Sep 21, 2011 at 4:47 PM, Tijmen Klein <T.R.Klein at student.rug.nl>
> wrote:****
>
> Hi everyone,****
>
> ** **
>
> I'm trying to get the end locations of the seeds inserted into
> a vtkStreamTracer. I insert some particles using a vtkDataSet, and I want to
> know the locations of where these particles ended up. How is this possible
> with a vtkStreamTracer?****
>
> ** **
>
> Some context: I'm working with time dependent data, and I want to use the
> end locations as the new seed locations for the next timestep. I only want
> to visualize the path (using streamtubes) of the current timestep.****
>
> ** **
>
> Cheers,****
>
> Tijmen****
>
> ** **
>
> _______________________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the VTK FAQ at:
> http://www.vtk.org/Wiki/VTK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.vtk.org/mailman/listinfo/vtkusers
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20111003/a134061b/attachment.htm>
More information about the vtkusers
mailing list