1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkGeoProjectionSource.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 "vtkGeoProjectionSource.h"
21
22 #include "vtkActor.h"
23 #include "vtkBox.h"
24 #include "vtkCellArray.h"
25 #include "vtkClipPolyData.h"
26 #include "vtkDoubleArray.h"
27 #include "vtkGeoGraticule.h"
28 #include "vtkGeoProjection.h"
29 #include "vtkGeoTerrainNode.h"
30 #include "vtkGeoTransform.h"
31 #include "vtkImageData.h"
32 #include "vtkMath.h"
33 #include "vtkMutexLock.h"
34 #include "vtkObjectFactory.h"
35 #include "vtkPlane.h"
36 #include "vtkPlanes.h"
37 #include "vtkPointData.h"
38 #include "vtkPolyData.h"
39 #include "vtkPolyDataMapper.h"
40 #include "vtkPropCollection.h"
41 #include "vtkProperty.h"
42 #include "vtkRenderer.h"
43 #include "vtkSmartPointer.h"
44 #include "vtkTransform.h"
45 #include "vtkTransformFilter.h"
46
47 #include <stack>
48 #include <utility>
49
50
51 // Heap dump method used for debugging threading issues on windows.
52 //void heapdump( void )
53 // {
54 // _HEAPINFO hinfo;
55 // int heapstatus;
56 // hinfo._pentry = nullptr;
57 // while( ( heapstatus = _heapwalk( &hinfo ) ) == _HEAPOK )
58 // {
59 // //printf( "%6s block at %Fp of size %4.4X\n", ( hinfo._useflag == _USEDENTRY ? "USED" : "FREE" ), hinfo._pentry, hinfo._size );
60 // }
61 //
62 // switch( heapstatus )
63 // {
64 // case _HEAPEMPTY:
65 // //printf( "OK - empty heap\n" );
66 // break;
67 // case _HEAPEND:
68 // //printf( "OK - end of heap\n" );
69 // break;
70 // case _HEAPBADPTR:
71 // printf( "ERROR - bad pointer to heap\n" );
72 // //throw("asdf");
73 // break;
74 // case _HEAPBADBEGIN:
75 // printf( "ERROR - bad start of heap\n" );
76 // //throw("asdf");
77 // break;
78 // case _HEAPBADNODE:
79 // printf( "ERROR - bad node in heap\n" );
80 // //throw("asdf");
81 // break;
82 // }
83 // }
84
85
86 vtkStandardNewMacro(vtkGeoProjectionSource);
87 vtkCxxSetObjectMacro(vtkGeoProjectionSource, Transform, vtkAbstractTransform);
88 //----------------------------------------------------------------------------
vtkGeoProjectionSource()89 vtkGeoProjectionSource::vtkGeoProjectionSource()
90 {
91 VTK_LEGACY_BODY(vtkGeoProjectionSource::vtkGeoProjectionSource, "VTK 8.2");
92 this->Projection=0;
93 this->Transform = nullptr;
94 this->MinCellsPerNode = 20;
95
96 this->TransformLock = vtkMutexLock::New();
97 }
98
99 //----------------------------------------------------------------------------
~vtkGeoProjectionSource()100 vtkGeoProjectionSource::~vtkGeoProjectionSource()
101 {
102 this->TransformLock->Delete();
103 this->SetTransform(nullptr);
104 }
105
PrintSelf(ostream & os,vtkIndent indent)106 void vtkGeoProjectionSource::PrintSelf( ostream& os, vtkIndent indent )
107 {
108 this->Superclass::PrintSelf( os, indent );
109 os << indent << "Projection: " << this->Projection << "\n";
110 os << indent << "Transform: " << this->Transform << "\n";
111 os << indent << "MinCellsPerNode: " << this->MinCellsPerNode << "\n";
112 }
113
114 //----------------------------------------------------------------------------
RefineAndComputeError(vtkGeoTerrainNode * node)115 void vtkGeoProjectionSource::RefineAndComputeError(vtkGeoTerrainNode* node)
116 {
117 double* latRange = node->GetLatitudeRange();
118 double* lonRange = node->GetLongitudeRange();
119
120 int level = node->GetGraticuleLevel();
121 double latDelta = vtkGeoGraticule::GetLatitudeDelta(level);
122 double lonDelta = vtkGeoGraticule::GetLongitudeDelta(level);
123 while ((latRange[1] - latRange[0])*(lonRange[1] - lonRange[0])/(latDelta*lonDelta) < this->MinCellsPerNode)
124 {
125 ++level;
126 latDelta = vtkGeoGraticule::GetLatitudeDelta(level);
127 lonDelta = vtkGeoGraticule::GetLongitudeDelta(level);
128 }
129
130 vtkSmartPointer<vtkGeoGraticule> grat = vtkSmartPointer<vtkGeoGraticule>::New();
131 vtkSmartPointer<vtkGeoGraticule> refinedGrat = vtkSmartPointer<vtkGeoGraticule>::New();
132 vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
133 vtkSmartPointer<vtkGeoTransform> trans = vtkSmartPointer<vtkGeoTransform>::New();
134 vtkSmartPointer<vtkGeoProjection> proj = vtkSmartPointer<vtkGeoProjection>::New();
135 proj->SetName(vtkGeoProjection::GetProjectionName(this->Projection));
136 trans->SetDestinationProjection(proj);
137 transformFilter->SetTransform(trans);
138 grat->SetGeometryType(vtkGeoGraticule::QUADRILATERALS);
139 grat->SetLatitudeBounds(latRange);
140 grat->SetLongitudeBounds(lonRange);
141 refinedGrat->SetGeometryType(vtkGeoGraticule::QUADRILATERALS);
142
143 vtkSmartPointer<vtkPolyData> geom = vtkSmartPointer<vtkPolyData>::New();
144 vtkSmartPointer<vtkPolyData> refined = vtkSmartPointer<vtkPolyData>::New();
145
146 do
147 {
148 grat->SetLatitudeLevel(level);
149 grat->SetLongitudeLevel(level);
150 transformFilter->SetInputConnection(grat->GetOutputPort());
151 transformFilter->Update();
152 geom->DeepCopy(transformFilter->GetOutput());
153 refinedGrat->SetLatitudeLevel(level+1);
154 refinedGrat->SetLongitudeLevel(level+1);
155 double* curLatRange = geom->GetPointData()->GetArray("LatLong")->GetRange(0);
156 refinedGrat->SetLatitudeBounds(curLatRange);
157 double* curLonRange = geom->GetPointData()->GetArray("LatLong")->GetRange(1);
158 refinedGrat->SetLongitudeBounds(curLonRange);
159 transformFilter->SetInputConnection(refinedGrat->GetOutputPort());
160 transformFilter->Update();
161 refined->DeepCopy(transformFilter->GetOutput());
162 ++level;
163 } while (geom->GetNumberOfCells() < this->MinCellsPerNode &&
164 level < vtkGeoGraticule::NUMBER_OF_LEVELS);
165
166 node->SetGraticuleLevel(level);
167
168 // Compute grid size
169 vtkDataArray* latLonArr = geom->GetPointData()->GetArray("LatLong");
170 double firstLon = latLonArr->GetComponent(0, 1);
171 vtkIdType gridSize[2] = {1, 0};
172 while (latLonArr->GetComponent(gridSize[0], 1) != firstLon)
173 {
174 gridSize[0]++;
175 }
176 gridSize[1] = geom->GetNumberOfPoints() / gridSize[0];
177
178 // Compute refined grid size
179 latLonArr = refined->GetPointData()->GetArray("LatLong");
180 firstLon = latLonArr->GetComponent(0, 1);
181 vtkIdType rgridSize[2] = {1, 0};
182 while (latLonArr->GetComponent(rgridSize[0], 1) != firstLon)
183 {
184 rgridSize[0]++;
185 }
186 rgridSize[1] = refined->GetNumberOfPoints() / rgridSize[0];
187
188 // Calculate error.
189 double error = 0.0;
190 double pt00[3];
191 double pt01[3];
192 double pt11[3];
193 double pt10[3];
194 double latFrac;
195 double lonFrac;
196 double curPt[3];
197 double interpPt[3];
198 double interpLon0;
199 double interpLon1;
200 double curError;
201
202 vtkIdType skip = (rgridSize[0]-1) / (gridSize[0]-1);
203 for (vtkIdType latInd = 0; latInd < rgridSize[1]-skip; ++latInd)
204 {
205 for (vtkIdType lonInd = 0; lonInd < rgridSize[0]-skip; ++lonInd)
206 {
207 vtkIdType ind00 = latInd * rgridSize[0] + lonInd;
208 vtkIdType ind01 = latInd * rgridSize[0] + lonInd + skip;
209 vtkIdType ind11 = ( latInd + skip ) * rgridSize[0] + lonInd + skip;
210 vtkIdType ind10 = ( latInd + skip ) * rgridSize[0] + lonInd;
211
212 refined->GetPoint(ind00, pt00);
213 refined->GetPoint(ind01, pt01);
214 refined->GetPoint(ind11, pt11);
215 refined->GetPoint(ind10, pt10);
216 for (vtkIdType rlatInd = latInd + 1; rlatInd < latInd + skip; ++rlatInd)
217 {
218 latFrac = static_cast<double>(rlatInd - latInd) / skip;
219 for (vtkIdType rlonInd = lonInd + 1; rlonInd < lonInd + skip; ++rlonInd)
220 {
221 lonFrac = static_cast<double>(rlonInd - lonInd) / skip;
222 refined->GetPoint(rlatInd * rgridSize[0] + rlonInd, curPt);
223 for (int c = 0; c < 3; ++c)
224 {
225 interpLon0 = (1.0 - lonFrac)*pt00[c] + lonFrac*pt01[c];
226 interpLon1 = (1.0 - lonFrac)*pt10[c] + lonFrac*pt11[c];
227 interpPt[c] = (1.0 - latFrac)*interpLon0 + latFrac*interpLon1;
228 }
229 curError = vtkMath::Distance2BetweenPoints(curPt, interpPt);
230 if (curError > error)
231 {
232 error = curError;
233 }
234 }
235 }
236 }
237 }
238
239 node->GetModel()->DeepCopy(geom);
240 node->SetError(sqrt(error));
241 }
242
243 //----------------------------------------------------------------------------
FetchRoot(vtkGeoTreeNode * r)244 bool vtkGeoProjectionSource::FetchRoot(vtkGeoTreeNode* r)
245 {
246 this->TransformLock->Lock();
247
248 vtkGeoTerrainNode* root = nullptr;
249 if (!(root = vtkGeoTerrainNode::SafeDownCast(r)))
250 {
251 vtkErrorMacro(<< "Can only fetch surface nodes from this source.");
252 return false;
253 }
254
255 // Let's start with graticule level 2 ... why not?
256 root->SetGraticuleLevel(2);
257
258 vtkSmartPointer<vtkGeoGraticule> grat = vtkSmartPointer<vtkGeoGraticule>::New();
259 grat->SetLatitudeLevel(root->GetGraticuleLevel());
260 grat->SetLongitudeLevel(root->GetGraticuleLevel());
261 grat->SetLongitudeBounds(-180.0, 180.0);
262 grat->SetLatitudeBounds(-90.0, 90.0);
263 grat->SetGeometryType(vtkGeoGraticule::QUADRILATERALS);
264
265 vtkSmartPointer<vtkTransformFilter> transformFilter = vtkSmartPointer<vtkTransformFilter>::New();
266 vtkSmartPointer<vtkGeoTransform> trans = vtkSmartPointer<vtkGeoTransform>::New();
267 vtkSmartPointer<vtkGeoProjection> proj = vtkSmartPointer<vtkGeoProjection>::New();
268 proj->SetName(vtkGeoProjection::GetProjectionName(this->Projection));
269 trans->SetDestinationProjection(proj);
270 transformFilter->SetTransform(trans);
271
272 transformFilter->SetInputConnection(grat->GetOutputPort());
273 transformFilter->Update();
274 const double *realBounds = transformFilter->GetOutput()->GetBounds();
275
276 // Extend the bounds just a tad to be safe
277 double bounds[4];
278 bounds[0] = realBounds[0] - (realBounds[1] - realBounds[0])*0.01;
279 bounds[1] = realBounds[1] + (realBounds[1] - realBounds[0])*0.01;
280 bounds[2] = realBounds[2] - (realBounds[3] - realBounds[2])*0.01;
281 bounds[3] = realBounds[3] + (realBounds[3] - realBounds[2])*0.01;
282
283 // Make the bounds square
284 if (bounds[1] - bounds[0] > bounds[3] - bounds[2])
285 {
286 double size = bounds[1] - bounds[0];
287 double center = (bounds[2] + bounds[3])/2.0;
288 bounds[2] = center - size/2.0;
289 bounds[3] = center + size/2.0;
290 }
291 else
292 {
293 double size = bounds[3] - bounds[2];
294 double center = (bounds[0] + bounds[1])/2.0;
295 bounds[0] = center - size/2.0;
296 bounds[1] = center + size/2.0;
297 }
298
299 root->GetModel()->DeepCopy(transformFilter->GetOutput());
300 root->SetLatitudeRange(-90.0, 90.0);
301 root->SetLongitudeRange(-180.0, 180.0);
302 root->SetProjectionBounds(bounds);
303 root->SetLevel(0);
304 this->RefineAndComputeError(root);
305
306 // Make sure bounds are up to date so we don't have threading issues
307 // when we hand this off to the main thread.
308 root->GetModel()->ComputeBounds();
309
310 this->TransformLock->Unlock();
311 return true;
312 }
313
314 //----------------------------------------------------------------------------
FetchChild(vtkGeoTreeNode * p,int index,vtkGeoTreeNode * c)315 bool vtkGeoProjectionSource::FetchChild(vtkGeoTreeNode* p, int index, vtkGeoTreeNode* c)
316 {
317 this->TransformLock->Lock();
318
319 vtkGeoTerrainNode* parent = nullptr;
320 if (!(parent = vtkGeoTerrainNode::SafeDownCast(p)))
321 {
322 vtkErrorMacro(<< "Can only fetch surface nodes from this source.");
323 return false;
324 }
325 vtkGeoTerrainNode* child = nullptr;
326 if (!(child = vtkGeoTerrainNode::SafeDownCast(c)))
327 {
328 vtkErrorMacro(<< "Can only fetch surface nodes from this source.");
329 return false;
330 }
331 if (!parent->HasData())
332 {
333 return false;
334 }
335
336 // Clip the cells of the children
337 double bounds[4];
338 parent->GetProjectionBounds(bounds);
339 double center[3] = {
340 (bounds[1] + bounds[0])/2.0,
341 (bounds[3] + bounds[2])/2.0, 0.0};
342
343 vtkSmartPointer<vtkClipPolyData> lonClip = vtkSmartPointer<vtkClipPolyData>::New();
344 vtkSmartPointer<vtkPlane> lonClipPlane = vtkSmartPointer<vtkPlane>::New();
345 lonClipPlane->SetOrigin(center);
346 lonClipPlane->SetNormal(-1.0, 0.0, 0.0);
347 lonClip->SetClipFunction(lonClipPlane);
348 lonClip->GenerateClippedOutputOn();
349 lonClip->SetInputData(parent->GetModel());
350
351 vtkSmartPointer<vtkPlane> latClipPlane = vtkSmartPointer<vtkPlane>::New();
352 latClipPlane->SetOrigin(center);
353 latClipPlane->SetNormal(0.0, -1.0, 0.0);
354 vtkSmartPointer<vtkClipPolyData> latClip = vtkSmartPointer<vtkClipPolyData>::New();
355 latClip->SetClipFunction(latClipPlane);
356 latClip->GenerateClippedOutputOn();
357 if (index % 2)
358 {
359 latClip->SetInputConnection(lonClip->GetOutputPort(1));
360 bounds[0] = center[0];
361 }
362 else
363 {
364 latClip->SetInputConnection(lonClip->GetOutputPort(0));
365 bounds[1] = center[0];
366 }
367 latClip->Update();
368 if (index / 2)
369 {
370 child->GetModel()->DeepCopy(latClip->GetOutput(1));
371 bounds[2] = center[1];
372 }
373 else
374 {
375 child->GetModel()->DeepCopy(latClip->GetOutput(0));
376 bounds[3] = center[1];
377 }
378 int level = parent->GetLevel() + 1;
379 child->SetLevel(level);
380 child->SetProjectionBounds(bounds);
381
382 // Set the id
383 if (level <= 15)
384 {
385 int id = parent->GetId() | (index << (level*2 - 2));
386 child->SetId(id);
387 }
388
389 double* latRange = nullptr;
390 double* lonRange = nullptr;
391 if (child->GetModel()->GetNumberOfPoints() > 0)
392 {
393 latRange = child->GetModel()->GetPointData()->GetArray("LatLong")->GetRange(0);
394 latRange[0] = (latRange[0] < -90) ? -90 : latRange[0];
395 latRange[1] = (latRange[1] > 90) ? 90 : latRange[1];
396 child->SetLatitudeRange(latRange);
397 lonRange = child->GetModel()->GetPointData()->GetArray("LatLong")->GetRange(1);
398 lonRange[0] = (lonRange[0] < -180) ? -180 : lonRange[0];
399 lonRange[1] = (lonRange[1] > 180) ? 180 : lonRange[1];
400 child->SetLongitudeRange(lonRange);
401 }
402 else
403 {
404 child->SetLatitudeRange(0.0, 0.0);
405 child->SetLongitudeRange(0.0, 0.0);
406 this->TransformLock->Unlock();
407 return true;
408 }
409
410 // Start with at least graticule level 2.
411 child->SetGraticuleLevel(2);
412
413 // Refine the node using vtkGeoGraticule and compute the error of the node.
414 this->RefineAndComputeError(child);
415
416 // We need to do four planar clips to get the desired result.
417 // Using vtkBox or vtkPlanes produces a fuzzy clip that is not acceptable.
418 for (int i = 0; i < 4; ++i)
419 {
420 vtkSmartPointer<vtkClipPolyData> finalClip = vtkSmartPointer<vtkClipPolyData>::New();
421 vtkSmartPointer<vtkPlane> plane = vtkSmartPointer<vtkPlane>::New();
422 if (i == 0)
423 {
424 plane->SetOrigin(bounds[0], 0.0, 0.0);
425 plane->SetNormal(1.0, 0.0, 0.0);
426 }
427 else if (i == 1)
428 {
429 plane->SetOrigin(bounds[1], 0.0, 0.0);
430 plane->SetNormal(-1.0, 0.0, 0.0);
431 }
432 else if (i == 2)
433 {
434 plane->SetOrigin(0.0, bounds[2], 0.0);
435 plane->SetNormal(0.0, 1.0, 0.0);
436 }
437 else if (i == 3)
438 {
439 plane->SetOrigin(0.0, bounds[3], 0.0);
440 plane->SetNormal(0.0, -1.0, 0.0);
441 }
442 finalClip->SetClipFunction(plane);
443 vtkSmartPointer<vtkPolyData> pd = vtkSmartPointer<vtkPolyData>::New();
444 pd->DeepCopy(child->GetModel());
445 finalClip->SetInputData(pd);
446 finalClip->Update();
447 child->GetModel()->DeepCopy(finalClip->GetOutput());
448 }
449
450 // The lat/long range could have changed
451 if (child->GetModel()->GetNumberOfPoints() > 0)
452 {
453 latRange = child->GetModel()->GetPointData()->GetArray("LatLong")->GetRange(0);
454 latRange[0] = (latRange[0] < -90) ? -90 : latRange[0];
455 latRange[1] = (latRange[1] > 90) ? 90 : latRange[1];
456 child->SetLatitudeRange(latRange);
457 lonRange = child->GetModel()->GetPointData()->GetArray("LatLong")->GetRange(1);
458 lonRange[0] = (lonRange[0] < -180) ? -180 : lonRange[0];
459 lonRange[1] = (lonRange[1] > 180) ? 180 : lonRange[1];
460 child->SetLongitudeRange(lonRange);
461 }
462 else
463 {
464 child->SetLatitudeRange(0.0, 0.0);
465 child->SetLongitudeRange(0.0, 0.0);
466 }
467
468 // Make sure bounds are up to date so we don't have threading issues
469 // when we hand this off to the main thread.
470 child->GetModel()->ComputeBounds();
471
472 this->TransformLock->Unlock();
473 return true;
474 }
475
476 //----------------------------------------------------------------------------
SetProjection(int projection)477 void vtkGeoProjectionSource::SetProjection(int projection)
478 {
479 this->Projection = projection;
480 vtkSmartPointer<vtkGeoTransform> trans = vtkSmartPointer<vtkGeoTransform>::New();
481 vtkSmartPointer<vtkGeoProjection> proj = vtkSmartPointer<vtkGeoProjection>::New();
482 proj->SetName(vtkGeoProjection::GetProjectionName(projection));
483 trans->SetDestinationProjection(proj);
484 this->SetTransform(trans);
485 }
486
487 //----------------------------------------------------------------------------
GetTransform()488 vtkAbstractTransform* vtkGeoProjectionSource::GetTransform()
489 {
490 return this->Transform;
491 }
492