vtkGlobeSource.cxx 12 KB
Newer Older
Ken Martin's avatar
Ken Martin committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGlobeSource.cxx

  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.

=========================================================================*/

/*-------------------------------------------------------------------------
  Copyright 2008 Sandia Corporation.
  Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
  the U.S. Government retains certain rights in this software.
-------------------------------------------------------------------------*/

#include "vtkGlobeSource.h"

#include "vtkCellArray.h"
25
#include "vtkDoubleArray.h"
Ken Martin's avatar
Ken Martin committed
26 27 28 29 30 31 32 33 34 35 36 37
#include "vtkFloatArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkGeoMath.h"
#include "vtkTimerLog.h"

Sean McBride's avatar
Sean McBride committed
38
#include <cmath>
Ken Martin's avatar
Ken Martin committed
39 40 41 42 43 44 45 46

vtkStandardNewMacro(vtkGlobeSource);

  // 0=NE, 1=SE, 2=SW, 3=NW

//----------------------------------------------------------------------------
vtkGlobeSource::vtkGlobeSource()
{
Sankhesh Jhaveri's avatar
Sankhesh Jhaveri committed
47
  VTK_LEGACY_BODY(vtkGlobeSource::vtkGlobeSource, "VTK 8.2");
48
  this->Origin[0] = this->Origin[1] = this->Origin[2] = 0.0;
49
  this->Radius = vtkGeoMath::EarthRadiusMeters();
50
  this->AutoCalculateCurtainHeight = true;
Ken Martin's avatar
Ken Martin committed
51
  this->CurtainHeight = 1000.0;
52

Ken Martin's avatar
Ken Martin committed
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
  this->LongitudeResolution = 10;
  this->LatitudeResolution = 10;
  this->StartLongitude = 0.0;
  this->EndLongitude = 360.0;
  this->StartLatitude = 0.0;
  this->EndLatitude = 180.0;
  this->QuadrilateralTessellation = 0;

  this->SetNumberOfInputPorts(0);
}

//----------------------------------------------------------------------------
void vtkGlobeSource::ComputeGlobePoint(
  double theta, double phi, double radius, double* x, double* normal)
{
68 69 70 71 72
  // Lets keep this conversion code in a single place.
  double tmp = cos( vtkMath::RadiansFromDegrees( phi ) );
  double n0 = -tmp * sin( vtkMath::RadiansFromDegrees( theta ) );
  double n1 = tmp * cos( vtkMath::RadiansFromDegrees( theta ) );
  double n2 = sin( vtkMath::RadiansFromDegrees( phi ) );
73

74 75 76
  x[0] = n0 * radius;
  x[1] = n1 * radius;
  x[2] = n2 * radius;
77

78
  if (normal)
79
  {
80 81 82
    normal[0] = n0;
    normal[1] = n1;
    normal[2] = n2;
83
  }
Ken Martin's avatar
Ken Martin committed
84 85 86 87 88 89 90 91 92 93
}

//----------------------------------------------------------------------------
void vtkGlobeSource::ComputeLatitudeLongitude(
  double* x, double& theta, double& phi)
{
  double rho = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
  double S = sqrt(x[0]*x[0] + x[1]*x[1]);
  phi = acos(x[2] / rho);
  if (x[0] >= 0)
94
  {
Ken Martin's avatar
Ken Martin committed
95
    theta = asin(x[1] / S);
96
  }
Ken Martin's avatar
Ken Martin committed
97
  else
98
  {
Ken Martin's avatar
Ken Martin committed
99
    theta = vtkMath::Pi() - asin(x[1] / S);
100
  }
101 102
  phi =   vtkMath::DegreesFromRadians( vtkMath::Pi() / 2.0 - phi );
  theta = vtkMath::DegreesFromRadians( theta - vtkMath::Pi()/2.0 );
Ken Martin's avatar
Ken Martin committed
103 104 105 106 107 108
}

//----------------------------------------------------------------------------
void vtkGlobeSource::AddPoint(
  double theta, double phi, double radius,
  vtkPoints* newPoints, vtkFloatArray* newNormals,
109 110
  vtkFloatArray* newLongitudeArray, vtkFloatArray* newLatitudeArray,
  vtkDoubleArray* newLatLongArray)
Ken Martin's avatar
Ken Martin committed
111 112
{
  double x[3], n[3];
113

Ken Martin's avatar
Ken Martin committed
114
  vtkGlobeSource::ComputeGlobePoint(theta, phi, radius, x, n);
115

116 117 118
  x[0] -= this->Origin[0];
  x[1] -= this->Origin[1];
  x[2] -= this->Origin[2];
119

Ken Martin's avatar
Ken Martin committed
120 121
  newPoints->InsertNextPoint(x);
  newNormals->InsertNextTuple(n);
122

Ken Martin's avatar
Ken Martin committed
123 124
  newLongitudeArray->InsertNextValue(theta);
  newLatitudeArray->InsertNextValue(phi);
125 126
  newLatLongArray->InsertNextValue(phi);
  newLatLongArray->InsertNextValue(theta);
Ken Martin's avatar
Ken Martin committed
127 128 129 130 131 132 133 134 135 136 137 138 139
}



//----------------------------------------------------------------------------
int vtkGlobeSource::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *outputVector)
{
  // get the info object
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

140
  // I used this to see if the background thread was really
Ken Martin's avatar
Ken Martin committed
141 142 143 144 145
  // operating asynchronously.
  //Sleep(2000);

  // I am going to compute the curtain height based on the level of the
  // terrain patch.
146
  if(this->AutoCalculateCurtainHeight)
147
  {
148 149
    this->CurtainHeight = (this->EndLongitude-this->StartLongitude)
      * this->Radius / 3600.0;
150
  }
Ken Martin's avatar
Ken Martin committed
151

luz.paz's avatar
luz.paz committed
152
  // get the output
Ken Martin's avatar
Ken Martin committed
153 154 155 156 157 158 159 160
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

  int i, j;
  int numPts, numPolys;
  vtkPoints *newPoints;
  vtkFloatArray *newLongitudeArray;
  vtkFloatArray *newLatitudeArray;
161
  vtkDoubleArray *newLatLongArray;
Ken Martin's avatar
Ken Martin committed
162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
  vtkFloatArray *newNormals;
  vtkCellArray *newPolys;
  double phi, theta;
  vtkIdType pts[4];

  // Set things up; allocate memory
  //
  numPts = this->LatitudeResolution * this->LongitudeResolution;
  // creating triangles
  numPolys = (this->LatitudeResolution-1)*(this->LongitudeResolution-1) * 2;

  // Add more for curtains.
  numPts += 2*(this->LatitudeResolution + this->LongitudeResolution);

  newPoints = vtkPoints::New();
  newPoints->Allocate(numPts);
  newNormals = vtkFloatArray::New();
  newNormals->SetNumberOfComponents(3);
  newNormals->Allocate(3*numPts);
  newNormals->SetName("Normals");

  newLongitudeArray = vtkFloatArray::New();
  newLongitudeArray->SetNumberOfComponents(1);
  newLongitudeArray->Allocate(numPts);
  newLongitudeArray->SetName("Longitude");

  newLatitudeArray = vtkFloatArray::New();
  newLatitudeArray->SetNumberOfComponents(1);
  newLatitudeArray->Allocate(numPts);
  newLatitudeArray->SetName("Latitude");
192 193 194 195 196

  newLatLongArray = vtkDoubleArray::New();
  newLatLongArray->SetNumberOfComponents(2);
  newLatLongArray->Allocate(2*numPts);
  newLatLongArray->SetName("LatLong");
197

Ken Martin's avatar
Ken Martin committed
198 199 200 201 202 203 204 205 206 207
  newPolys = vtkCellArray::New();
  newPolys->Allocate(newPolys->EstimateSize(numPolys, 3));

  // Create sphere
  //

  double deltaLongitude;
  double deltaLatitude;

  // Check data, determine increments, and convert to radians
208
  deltaLongitude = (this->EndLongitude - this->StartLongitude)
209
                 / static_cast<double>(this->LongitudeResolution-1);
Ken Martin's avatar
Ken Martin committed
210

211
  deltaLatitude = (this->EndLatitude - this->StartLatitude)
212
                 / static_cast<double>(this->LatitudeResolution-1);
213

Ken Martin's avatar
Ken Martin committed
214 215
  // Create points and point data.
  for (j=0; j<this->LatitudeResolution; j++)
216
  {
217
    phi = this->StartLatitude + j*deltaLatitude;
Ken Martin's avatar
Ken Martin committed
218
    for (i=0; i < this->LongitudeResolution; i++)
219
    {
220
      theta = this->StartLongitude + i*deltaLongitude;
221
      this->AddPoint(theta, phi, this->Radius,
Ken Martin's avatar
Ken Martin committed
222
                     newPoints, newNormals,
223 224
                     newLongitudeArray, newLatitudeArray,
                     newLatLongArray);
225
    }
Ken Martin's avatar
Ken Martin committed
226 227
    this->UpdateProgress(
      0.10 + 0.50*j/static_cast<float>(this->LatitudeResolution));
228
  }
Ken Martin's avatar
Ken Martin committed
229 230 231

  // Create the extra points for the curtains.
  for (i=0; i < this->LongitudeResolution; i++)
232
  {
233
    theta = this->StartLongitude + i*deltaLongitude;
Ken Martin's avatar
Ken Martin committed
234
    phi = this->StartLatitude;
235
    this->AddPoint(theta, phi, this->Radius-this->CurtainHeight,
Ken Martin's avatar
Ken Martin committed
236
                   newPoints, newNormals,
237 238
                   newLongitudeArray, newLatitudeArray,
                   newLatLongArray);
239
  }
Ken Martin's avatar
Ken Martin committed
240
  for (i=0; i < this->LongitudeResolution; i++)
241
  {
242
    theta = this->StartLongitude + i*deltaLongitude;
Ken Martin's avatar
Ken Martin committed
243
    phi = this->EndLatitude;
244
    this->AddPoint(theta, phi, this->Radius-this->CurtainHeight,
Ken Martin's avatar
Ken Martin committed
245
                   newPoints, newNormals,
246
                   newLongitudeArray, newLatitudeArray,
247
                   newLatLongArray);
248
  }
Ken Martin's avatar
Ken Martin committed
249
  for (j=0; j < this->LatitudeResolution; j++)
250
  {
Ken Martin's avatar
Ken Martin committed
251
    theta = this->StartLongitude;
252
    phi = this->StartLatitude + j*deltaLatitude;
253
    this->AddPoint(theta, phi, this->Radius-this->CurtainHeight,
Ken Martin's avatar
Ken Martin committed
254
                   newPoints, newNormals,
255 256
                   newLongitudeArray, newLatitudeArray,
                   newLatLongArray);
257
  }
Ken Martin's avatar
Ken Martin committed
258
  for (j=0; j < this->LatitudeResolution; j++)
259
  {
Ken Martin's avatar
Ken Martin committed
260
    theta = this->EndLongitude;
261
    phi = this->StartLatitude + j*deltaLatitude;
262
    this->AddPoint(theta, phi, this->Radius-this->CurtainHeight,
Ken Martin's avatar
Ken Martin committed
263
                   newPoints, newNormals,
264
                   newLongitudeArray, newLatitudeArray,
265
                   newLatLongArray);
266
  }
267 268


Ken Martin's avatar
Ken Martin committed
269 270 271 272
  // Generate mesh connectivity
  vtkIdType rowId = 0;
  vtkIdType cornerId;
  for (j=1; j < this->LatitudeResolution; ++j)
273
  {
Ken Martin's avatar
Ken Martin committed
274 275
    cornerId = rowId;
    for (i=1; i < this->LongitudeResolution; ++i)
276
    {
Ken Martin's avatar
Ken Martin committed
277 278 279 280 281 282 283 284
      pts[0] = cornerId;
      pts[2] = cornerId + this->LongitudeResolution;
      pts[1] = pts[2] + 1;
      newPolys->InsertNextCell(3, pts);
      pts[2] = pts[1];
      pts[1] = cornerId + 1;
      newPolys->InsertNextCell(3, pts);
      ++cornerId;
285
    }
Ken Martin's avatar
Ken Martin committed
286 287 288
    rowId += this->LongitudeResolution;
    this->UpdateProgress(
      0.70 + 0.30*j/static_cast<double>(this->LatitudeResolution));
289
  }
290

Ken Martin's avatar
Ken Martin committed
291
  // Create curtain quads.
292
  vtkIdType curtainPointId =
Ken Martin's avatar
Ken Martin committed
293 294 295 296
    this->LongitudeResolution * this->LatitudeResolution;
  vtkIdType edgeOffset;
  edgeOffset = 0;
  for (i=1; i < this->LongitudeResolution; ++i)
297
  {
Ken Martin's avatar
Ken Martin committed
298 299 300 301 302 303
    pts[0] = edgeOffset + i; // i starts at 1.
    pts[1] = pts[0] - 1;
    pts[2] = curtainPointId;
    pts[3] = curtainPointId + 1;
    newPolys->InsertNextCell(4, pts);
    ++curtainPointId;
304
  }
Ken Martin's avatar
Ken Martin committed
305 306 307
  ++curtainPointId; // Skip 2 to the next edge.
  edgeOffset = (this->LongitudeResolution)*(this->LatitudeResolution-1);
  for (i=1; i < this->LongitudeResolution; ++i)
308
  {
Ken Martin's avatar
Ken Martin committed
309 310 311 312 313 314
    pts[0] = edgeOffset + i - 1; // i starts at 1
    pts[1] = pts[0] + 1;
    pts[2] = curtainPointId + 1;
    pts[3] = curtainPointId;
    newPolys->InsertNextCell(4, pts);
    ++curtainPointId;
315
  }
Ken Martin's avatar
Ken Martin committed
316 317 318
  ++curtainPointId;
  edgeOffset = 0;
  for (j=1; j < this->LatitudeResolution; ++j)
319
  {
Ken Martin's avatar
Ken Martin committed
320 321 322 323 324 325
    pts[0] = edgeOffset + j*this->LongitudeResolution;
    pts[1] = pts[0] - this->LongitudeResolution;
    pts[2] = curtainPointId;
    pts[3] = curtainPointId + 1;
    newPolys->InsertNextCell(4, pts);
    ++curtainPointId;
326
  }
Ken Martin's avatar
Ken Martin committed
327 328 329
  ++curtainPointId;
  edgeOffset = (this->LongitudeResolution-1);
  for (j=1; j < this->LatitudeResolution; ++j)
330
  {
Ken Martin's avatar
Ken Martin committed
331 332 333 334 335 336
    pts[0] = edgeOffset + (j-1)*this->LongitudeResolution;
    pts[1] = pts[0] + this->LongitudeResolution;
    pts[2] = curtainPointId + 1;
    pts[3] = curtainPointId;
    newPolys->InsertNextCell(4, pts);
    ++curtainPointId;
337
  }
338 339


luz.paz's avatar
luz.paz committed
340
  // Update ourselves and release memory
Ken Martin's avatar
Ken Martin committed
341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357
  //
  newPoints->Squeeze();
  output->SetPoints(newPoints);
  newPoints->Delete();

  newNormals->Squeeze();
  output->GetPointData()->SetNormals(newNormals);
  newNormals->Delete();

  newLongitudeArray->Squeeze();
  output->GetPointData()->AddArray(newLongitudeArray);
  newLongitudeArray->Delete();

  newLatitudeArray->Squeeze();
  output->GetPointData()->AddArray(newLatitudeArray);
  newLatitudeArray->Delete();

358 359 360 361
  newLatLongArray->Squeeze();
  output->GetPointData()->AddArray(newLatLongArray);
  newLatLongArray->Delete();

Ken Martin's avatar
Ken Martin committed
362 363 364 365 366 367 368 369 370 371 372 373
  newPolys->Squeeze();
  output->SetPolys(newPolys);
  newPolys->Delete();

  return 1;
}

//----------------------------------------------------------------------------
void vtkGlobeSource::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);

374 375
  os << indent << "AutoCalculateCurtainHeight: " <<
    (this->AutoCalculateCurtainHeight ? "ON" : "OFF") << "\n";
Ken Martin's avatar
Ken Martin committed
376
  os << indent << "CurtainHeight: " << this->CurtainHeight << "\n";
Ken Martin's avatar
Ken Martin committed
377 378 379 380 381 382 383
  os << indent << "Longitude Resolution: " << this->LongitudeResolution << "\n";
  os << indent << "Latitude Resolution: " << this->LatitudeResolution << "\n";
  os << indent << "Longitude Start: " << this->StartLongitude << "\n";
  os << indent << "Latitude Start: " << this->StartLatitude << "\n";
  os << indent << "Longitude End: " << this->EndLongitude << "\n";
  os << indent << "Latitude End: " << this->EndLatitude << "\n";
  os << indent << "Radius: " << this->Radius << "\n";
384 385 386
  os << indent << "Origin: " << this->Origin[0] << ","
                             << this->Origin[1] << ","
                             << this->Origin[2] << "\n";
387 388
  os << indent
     << "Quadrilateral Tessellation: "
Ken Martin's avatar
Ken Martin committed
389 390
     << this->QuadrilateralTessellation << "\n";
}