1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkStreamLine.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 "vtkStreamLine.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkDataSet.h"
19 #include "vtkFloatArray.h"
20 #include "vtkFloatArray.h"
21 #include "vtkInformation.h"
22 #include "vtkInformationVector.h"
23 #include "vtkMath.h"
24 #include "vtkObjectFactory.h"
25 #include "vtkPointData.h"
26 #include "vtkPolyData.h"
27 #include "vtkPolyLine.h"
28
29 vtkStandardNewMacro(vtkStreamLine);
30
31 // Construct object with step size set to 1.0.
vtkStreamLine()32 vtkStreamLine::vtkStreamLine()
33 {
34 this->StepLength = 1.0;
35 this->NumberOfStreamers = 0;
36 }
37
RequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)38 int vtkStreamLine::RequestData(
39 vtkInformation *,
40 vtkInformationVector **inputVector,
41 vtkInformationVector *outputVector)
42 {
43 vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
44 vtkInformation *outInfo = outputVector->GetInformationObject(0);
45 vtkInformation *sourceInfo = inputVector[1]->GetInformationObject(0);
46
47 vtkDataSet *input = vtkDataSet::SafeDownCast(
48 inInfo->Get(vtkDataObject::DATA_OBJECT()));
49 vtkPolyData *output = vtkPolyData::SafeDownCast(
50 outInfo->Get(vtkDataObject::DATA_OBJECT()));
51 vtkDataSet *source = 0;
52 if (sourceInfo)
53 {
54 source = vtkDataSet::SafeDownCast(
55 sourceInfo->Get(vtkDataObject::DATA_OBJECT()));
56 }
57
58 vtkStreamer::StreamPoint *sPrev, *sPtr;
59 vtkPoints *newPts;
60 vtkFloatArray *newVectors;
61 vtkFloatArray *newScalars=NULL;
62 vtkCellArray *newLines;
63 vtkIdType ptId, i, id;
64 int j;
65 vtkIdList *pts;
66 double tOffset, x[3], v[3], s, r;
67 double theta;
68 vtkPolyLine* lineNormalGenerator = NULL;
69 vtkFloatArray* normals = NULL;
70 vtkFloatArray* rotation = 0;
71
72 this->SavePointInterval = this->StepLength;
73 this->vtkStreamer::Integrate(input, source);
74 if ( this->NumberOfStreamers <= 0 ) {return 1;}
75
76 pts = vtkIdList::New();
77 pts->Allocate(2500);
78
79 //
80 // Convert streamer into lines. Lines may be dashed.
81 //
82 newPts = vtkPoints::New();
83 newPts->Allocate(1000);
84 newVectors = vtkFloatArray::New();
85 newVectors->SetNumberOfComponents(3);
86 newVectors->Allocate(3000);
87 if ( this->Vorticity )
88 {
89 lineNormalGenerator = vtkPolyLine::New();
90 normals = vtkFloatArray::New();
91 normals->SetNumberOfComponents(3);
92 normals->Allocate(3000);
93 rotation = vtkFloatArray::New();
94 rotation->SetNumberOfComponents(1);
95 rotation->Allocate(1000);
96 rotation->SetName("Thetas");
97 output->GetPointData()->AddArray(rotation);
98 }
99
100 if ( input->GetPointData()->GetScalars() || this->SpeedScalars
101 || this->OrientationScalars)
102 {
103 newScalars = vtkFloatArray::New();
104 newScalars->Allocate(1000);
105 }
106 newLines = vtkCellArray::New();
107 newLines->Allocate(newLines->EstimateSize(2*this->NumberOfStreamers,
108 VTK_CELL_SIZE));
109 //
110 // Loop over all streamers generating points
111 //
112 for (ptId=0; ptId < this->NumberOfStreamers; ptId++)
113 {
114 if ( this->Streamers[ptId].GetNumberOfPoints() < 2 )
115 {
116 continue;
117 }
118 sPrev = this->Streamers[ptId].GetStreamPoint(0);
119 sPtr = this->Streamers[ptId].GetStreamPoint(1);
120
121 if ( this->Streamers[ptId].GetNumberOfPoints() == 2 && sPtr->cellId >= 0 )
122 {
123 continue;
124 }
125
126 tOffset = sPrev->t;
127
128 for ( i=1;
129 i < this->Streamers[ptId].GetNumberOfPoints() && sPtr->cellId >= 0;
130 i++, sPrev=sPtr, sPtr=this->Streamers[ptId].GetStreamPoint(i) )
131 {
132 //
133 // Create points for line
134 //
135 while ( tOffset >= sPrev->t && tOffset < sPtr->t )
136 {
137 r = (tOffset - sPrev->t) / (sPtr->t - sPrev->t);
138
139 for (j=0; j<3; j++)
140 {
141 x[j] = sPrev->x[j] + r * (sPtr->x[j] - sPrev->x[j]);
142 v[j] = sPrev->v[j] + r * (sPtr->v[j] - sPrev->v[j]);
143 }
144
145 // add point to line
146 id = newPts->InsertNextPoint(x);
147 pts->InsertNextId(id);
148 newVectors->InsertTuple(id,v);
149
150 if ( newScalars )
151 {
152 s = sPrev->s + r * (sPtr->s - sPrev->s);
153 newScalars->InsertTuple(id,&s);
154 }
155
156 if ( this->Vorticity )
157 {
158 // Store the rotation values. Used after all the streamlines
159 // are generated.
160 theta = sPrev->theta + r * (sPtr->theta - sPrev->theta);
161 rotation->InsertTuple(id, &theta);
162 }
163
164 tOffset += this->StepLength;
165
166 } // while
167 } //for this streamer
168
169 if ( pts->GetNumberOfIds() > 1 )
170 {
171 newLines->InsertNextCell(pts);
172 pts->Reset();
173 }
174 } //for all streamers
175
176 vtkDebugMacro(<<"Created " << newPts->GetNumberOfPoints() << " points, "
177 << newLines->GetNumberOfCells() << " lines");
178
179 if (this->Vorticity)
180 {
181 // Rotate the normal vectors with stream vorticity
182 vtkIdType nPts=0;
183 vtkIdType *linePts=0;
184 double normal[3], local1[3], local2[3], length, costheta, sintheta;
185
186 lineNormalGenerator->GenerateSlidingNormals(newPts,newLines,normals);
187
188 // Loop over all lines, from the above code we are know that each line
189 // will have at least two points and that no points will be shared
190 // between lines. It is important to loop over the points used by the
191 // lines because newPts may actually contain points that are not used by
192 // any lines. The normals are only calculated for points that are used
193 // in lines so referencing normals for all points can lead to UMRs
194 for (newLines->InitTraversal(); newLines->GetNextCell(nPts,linePts); )
195 {
196 for(i=0; i<nPts; i++)
197 {
198 normals->GetTuple(linePts[i], normal);
199 newVectors->GetTuple(linePts[i], v);
200 // obtain two unit orthogonal vectors on the plane perpendicular to
201 // the streamline
202 for(j=0; j<3; j++)
203 {
204 local1[j] = normal[j];
205 }
206 length = vtkMath::Normalize(local1);
207 vtkMath::Cross(local1, v, local2);
208 vtkMath::Normalize(local2);
209 // Rotate the normal with theta
210 rotation->GetTuple(linePts[i], &theta);
211 costheta = cos(theta);
212 sintheta = sin(theta);
213 for(j=0; j<3; j++)
214 {
215 normal[j] = length* (costheta*local1[j] + sintheta*local2[j]);
216 }
217 normals->SetTuple(linePts[i], normal);
218 }
219 }
220 output->GetPointData()->SetNormals(normals);
221 normals->Delete();
222 lineNormalGenerator->Delete();
223 rotation->Delete();
224 }
225
226 output->SetPoints(newPts);
227 newPts->Delete();
228
229 output->GetPointData()->SetVectors(newVectors);
230 newVectors->Delete();
231
232 if ( newScalars )
233 {
234 int idx = output->GetPointData()->AddArray(newScalars);
235 output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
236 newScalars->Delete();
237 }
238
239 pts->Delete();
240 output->SetLines(newLines);
241 newLines->Delete();
242
243 // Delete the streamers since they are no longer needed
244 delete[] this->Streamers;
245 this->Streamers = 0;
246 this->NumberOfStreamers = 0;
247
248 output->Squeeze();
249
250 return 1;
251 }
252
PrintSelf(ostream & os,vtkIndent indent)253 void vtkStreamLine::PrintSelf(ostream& os, vtkIndent indent)
254 {
255 this->Superclass::PrintSelf(os,indent);
256
257 os << indent << "Step Length: " << this->StepLength << "\n";
258
259 }
260