1 /*=========================================================================
2 
3   Program:   Visualization Toolkit
4   Module:    vtkGeoAdaptiveArcs.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 /*-------------------------------------------------------------------------
16   Copyright 2008 Sandia Corporation.
17   Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
18   the U.S. Government retains certain rights in this software.
19 -------------------------------------------------------------------------*/
20 #include "vtkGeoAdaptiveArcs.h"
21 
22 #include "vtkCamera.h"
23 #include "vtkCellArray.h"
24 #include "vtkCellData.h"
25 #include "vtkCoordinate.h"
26 #include "vtkDoubleArray.h"
27 #include "vtkFloatArray.h"
28 #include "vtkGeoMath.h"
29 #include "vtkGlobeSource.h"
30 #include "vtkMath.h"
31 #include "vtkInformation.h"
32 #include "vtkInformationVector.h"
33 #include "vtkObjectFactory.h"
34 #include "vtkPointData.h"
35 #include "vtkRenderWindow.h"
36 #include "vtkRendererCollection.h"
37 #include "vtkTimerLog.h"
38 
39 #include <map>
40 
41 vtkStandardNewMacro(vtkGeoAdaptiveArcs);
42 
43 //-------------------------------------------------------------------------
vtkGeoAdaptiveArcs()44 vtkGeoAdaptiveArcs::vtkGeoAdaptiveArcs()
45 {
46   VTK_LEGACY_BODY(vtkGeoAdaptiveArcs::vtkGeoAdaptiveArcs, "VTK 8.2");
47   this->GlobeRadius = vtkGeoMath::EarthRadiusMeters();
48   this->Renderer = nullptr;
49   this->MaximumPixelSeparation = 10.0;
50   this->MinimumPixelSeparation = 1.0;
51   this->LastInput = nullptr;
52   this->LastInputMTime = 0;
53   this->InputLatitude = vtkDoubleArray::New();
54   this->InputLongitude = vtkDoubleArray::New();
55 }
56 
57 //-------------------------------------------------------------------------
~vtkGeoAdaptiveArcs()58 vtkGeoAdaptiveArcs::~vtkGeoAdaptiveArcs()
59 {
60   this->InputLatitude->Delete();
61   this->InputLongitude->Delete();
62 }
63 
64 //-------------------------------------------------------------------------
RequestData(vtkInformation * vtkNotUsed (request),vtkInformationVector ** inputVector,vtkInformationVector * outputVector)65 int vtkGeoAdaptiveArcs::RequestData(
66   vtkInformation *vtkNotUsed(request),
67   vtkInformationVector **inputVector,
68   vtkInformationVector *outputVector)
69 {
70   if (!this->Renderer)
71   {
72     vtkErrorMacro("Renderer cannot be null.");
73     return 0;
74   }
75   // get the info objects
76   vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
77   vtkInformation *outInfo = outputVector->GetInformationObject(0);
78 
79   // get the input and output
80   vtkPolyData *input = vtkPolyData::SafeDownCast(
81     inInfo->Get(vtkDataObject::DATA_OBJECT()));
82   vtkPolyData *output = vtkPolyData::SafeDownCast(
83     outInfo->Get(vtkDataObject::DATA_OBJECT()));
84 
85   // if the input has changed, compute helper arrays
86   if (input != this->LastInput ||
87       input->GetMTime() > this->LastInputMTime)
88   {
89     this->InputLatitude->Initialize();
90     this->InputLongitude->Initialize();
91     vtkPoints* points = input->GetPoints();
92     double curPtLL[2];
93     double curPoint[3];
94     for (vtkIdType i = 0; i < input->GetNumberOfPoints(); ++i)
95     {
96       points->GetPoint(i, curPoint);
97       vtkGlobeSource::ComputeLatitudeLongitude(
98         curPoint, curPtLL[0], curPtLL[1]);
99       this->InputLongitude->InsertNextValue(curPtLL[0]);
100       this->InputLatitude->InsertNextValue(curPtLL[1]);
101     }
102     this->LastInput = input;
103     this->LastInputMTime = input->GetMTime();
104   }
105 
106   // Traverse input lines, adding a circle for each line segment.
107   int *renSize = this->Renderer->GetSize();
108   // Maximum distance from center of renderer.
109   double maxDist = (renSize[0] > renSize[1]) ? 1.1*renSize[0]/2.0 : 1.1*renSize[1]/2.0;
110   vtkCellArray* lines = input->GetLines();
111   vtkCellArray* newLines = vtkCellArray::New();
112   vtkPoints* points = input->GetPoints();
113   float* pointsPtr = static_cast<float*>(points->GetVoidPointer(0));
114   vtkPoints* newPoints = vtkPoints::New();
115   double viewAngle = this->Renderer->GetActiveCamera()->GetViewAngle();
116   double cameraPos[3];
117   this->Renderer->GetActiveCamera()->GetPosition(cameraPos);
118   double cameraDir[3];
119   this->Renderer->GetActiveCamera()->GetDirectionOfProjection(cameraDir);
120   lines->InitTraversal();
121   for (vtkIdType i = 0; i < lines->GetNumberOfCells(); i++)
122   {
123     vtkIdType npts=0; // to remove warning
124     vtkIdType* pts=nullptr; // to remove warning
125     lines->GetNextCell(npts, pts);
126 
127     bool lastPointOffScreen = false;
128     bool lastPointTooClose = false;
129 #if defined(VTK_AGGRESSIVE_ARCS)
130     bool lastPointOnOtherSide = false;
131 #endif
132     double curPoint[3];
133     double lastPtLL[2] = {0.0, 0.0};
134     double curPtLL[2];
135     double lastVec[3] = {0.0, 0.0, 0.0};
136     double curVec[3];
137     curPoint[0] = pointsPtr[3*pts[0]+0];
138     curPoint[1] = pointsPtr[3*pts[0]+1];
139     curPoint[2] = pointsPtr[3*pts[0]+2];
140     curPtLL[0] = this->InputLongitude->GetValue(pts[0]);
141     curPtLL[1] = this->InputLatitude->GetValue(pts[0]);
142     double curVecSize = 0;
143     for (int c = 0; c < 3; ++c)
144     {
145       curVec[c] = curPoint[c] - cameraPos[c];
146       curVecSize += curVec[c]*curVec[c];
147     }
148     curVecSize = sqrt(curVecSize);
149     curVec[0] /= curVecSize;
150     curVec[1] /= curVecSize;
151     curVec[2] /= curVecSize;
152 
153     for (vtkIdType p = 1; p < npts; ++p)
154     {
155       // Advance the point unless the last point was too close.
156       if (!lastPointTooClose)
157       {
158         for (int c = 0; c < 3; ++c)
159         {
160 //          lastPoint[c] = curPoint[c];
161           lastVec[c] = curVec[c];
162         }
163         lastPtLL[0] = curPtLL[0];
164         lastPtLL[1] = curPtLL[1];
165       }
166 #if defined(VTK_AGGRESSIVE_ARCS)
167       // If this code is uncommented, then
168       // Be aggressive ... skip several points if the last
169       // one was offscreen or on the other side of the globe.
170       if (lastPointOffScreen || lastPointOnOtherSide)
171       {
172         p += 5;
173         if (p >= npts)
174         {
175           p = npts - 1;
176         }
177       }
178 #endif
179       curPoint[0] = pointsPtr[3*pts[p]+0];
180       curPoint[1] = pointsPtr[3*pts[p]+1];
181       curPoint[2] = pointsPtr[3*pts[p]+2];
182       curPtLL[0] = this->InputLongitude->GetValue(pts[p]);
183       curPtLL[1] = this->InputLatitude->GetValue(pts[p]);
184       curVecSize = 0;
185       for (int c = 0; c < 3; ++c)
186       {
187         curVec[c] = curPoint[c] - cameraPos[c];
188         curVecSize += curVec[c]*curVec[c];
189       }
190       curVecSize = sqrt(curVecSize);
191       curVec[0] /= curVecSize;
192       curVec[1] /= curVecSize;
193       curVec[2] /= curVecSize;
194 
195       // Clear flags
196 #if defined(VTK_AGGRESSIVE_ARCS)
197       lastPointOnOtherSide = false;
198 #endif
199       lastPointOffScreen = false;
200       lastPointTooClose = false;
201 
202       // Don't draw lines off the current screen.
203       double distFromCenterApprox =
204         vtkMath::DegreesFromRadians( acos( curVec[0] * cameraDir[0] +
205                                            curVec[1] * cameraDir[1] +
206                                            curVec[2] * cameraDir[2] ) )
207         / viewAngle * renSize[1];
208       if (distFromCenterApprox > maxDist)
209       {
210         // If both last point and this point are offscreen, skip
211         // drawing the line
212         if (lastPointOffScreen == true)
213         {
214           continue;
215         }
216         lastPointOffScreen = true;
217       }
218 
219       // Don't draw lines on the other side of the world.
220       //if (vtkMath::Dot(curPoint, cameraPos) < 0)
221       if (curPoint[0]*cameraPos[0]+
222           curPoint[1]*cameraPos[1]+
223           curPoint[2]*cameraPos[2] < 0)
224       {
225 #if defined(VTK_AGGRESSIVE_ARCS)
226         lastPointOnOtherSide = true;
227 #endif
228         continue;
229       }
230 
231       double distApprox =
232         vtkMath::DegreesFromRadians( acos( lastVec[0] * curVec[0] +
233                                            lastVec[1] * curVec[1] +
234                                            lastVec[2] * curVec[2] ) )
235         / viewAngle * renSize[1];
236 
237       // If the points are too close, skip over it to the next point.
238       if (distApprox < this->MinimumPixelSeparation)
239       {
240         lastPointTooClose = true;
241         continue;
242       }
243 
244       // Calculate the number of subdivisions.
245       vtkIdType numDivisions = static_cast<vtkIdType>(distApprox / this->MaximumPixelSeparation + 0.5) + 1;
246       if (numDivisions < 2)
247       {
248         numDivisions = 2;
249       }
250 
251       // Create the new cell
252       newLines->InsertNextCell(numDivisions);
253 
254       for (vtkIdType s = 0; s < numDivisions; ++s)
255       {
256         // Interpolate in lat-long.
257         double interpPtLL[2];
258         double frac = static_cast<double>(s) / (numDivisions - 1);
259         for (int c = 0; c < 2; ++c)
260         {
261           interpPtLL[c] = frac*curPtLL[c] + (1.0 - frac)*lastPtLL[c];
262         }
263         // Convert lat-long to world;
264         double interpPt[3];
265         vtkGlobeSource::ComputeGlobePoint(interpPtLL[0], interpPtLL[1], this->GlobeRadius, interpPt);
266         vtkIdType newPt = newPoints->InsertNextPoint(interpPt);
267         newLines->InsertCellPoint(newPt);
268       }
269 
270     }
271   }
272 
273   // Send the data to output.
274   output->SetLines(newLines);
275   output->SetPoints(newPoints);
276 
277   // Clean up.
278   newLines->Delete();
279   newPoints->Delete();
280 
281   return 1;
282 }
283 
284 //-------------------------------------------------------------------------
SetRenderer(vtkRenderer * ren)285 void vtkGeoAdaptiveArcs::SetRenderer(vtkRenderer *ren)
286 {
287   // Do not reference count this, it will cause a loop.
288   this->Renderer = ren;
289 }
290 
291 //-------------------------------------------------------------------------
GetMTime()292 vtkMTimeType vtkGeoAdaptiveArcs::GetMTime()
293 {
294   vtkMTimeType retMTime = this->Superclass::GetMTime();
295   if ( this->Renderer )
296   {
297     vtkMTimeType tmpTime = this->Renderer->GetMTime();
298     if ( tmpTime > retMTime )
299     {
300       retMTime = tmpTime;
301     }
302     vtkCamera* cam = this->Renderer->GetActiveCamera();
303     if ( cam )
304     {
305       tmpTime = cam->GetMTime();
306       if ( tmpTime > retMTime )
307        retMTime = tmpTime;
308     }
309   }
310   return retMTime;
311 }
312 
313 //-------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)314 void vtkGeoAdaptiveArcs::PrintSelf(ostream& os, vtkIndent indent)
315 {
316   this->Superclass::PrintSelf(os,indent);
317   os << indent << "GlobeRadius: " << this->GlobeRadius << endl;
318   os << indent << "MinimumPixelSeparation: " << this->MinimumPixelSeparation << endl;
319   os << indent << "MaximumPixelSeparation: " << this->MaximumPixelSeparation << endl;
320   os << indent << "Renderer: " << (this->Renderer ? "" : "(null)") << endl;
321   if (this->Renderer)
322   {
323     this->Renderer->PrintSelf(os, indent.GetNextIndent());
324   }
325 }
326 
327