vtkGeoArcs.cxx 5.73 KB
Newer Older
Ken Martin's avatar
Ken Martin committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkGeoArcs.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.

=========================================================================*/
15 16 17 18 19
/*-------------------------------------------------------------------------
  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.
-------------------------------------------------------------------------*/
Ken Martin's avatar
Ken Martin committed
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
#include "vtkGeoArcs.h"

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkFloatArray.h"
#include "vtkGeoMath.h"
#include "vtkGlobeSource.h"
#include "vtkMath.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"

vtkStandardNewMacro(vtkGeoArcs);

vtkGeoArcs::vtkGeoArcs()
{
  this->GlobeRadius = vtkGeoMath::EarthRadiusMeters();
  this->ExplodeFactor = 0.2;
  this->NumberOfSubdivisions = 20;
}

int vtkGeoArcs::RequestData(
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

51
  // get the input and output
Ken Martin's avatar
Ken Martin committed
52 53 54 55
  vtkPolyData *input = vtkPolyData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData *output = vtkPolyData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));
56

Ken Martin's avatar
Ken Martin committed
57 58 59 60 61 62 63 64 65 66
  // Prepare to copy cell data
  output->GetCellData()->CopyAllocate(input->GetCellData());

  // Traverse input lines, adding a circle for each line segment.
  vtkCellArray* lines = input->GetLines();
  vtkCellArray* newLines = vtkCellArray::New();
  vtkPoints* newPoints = vtkPoints::New();
  newPoints->DeepCopy(input->GetPoints());
  lines->InitTraversal();
  for (vtkIdType i = 0; i < lines->GetNumberOfCells(); i++)
67
  {
68
      vtkIdType npts=0; // to remove warning
69
    vtkIdType* pts=nullptr; // to remove warning
Ken Martin's avatar
Ken Martin committed
70
    lines->GetNextCell(npts, pts);
71

Ken Martin's avatar
Ken Martin committed
72 73
    double lastPoint[3];
    newPoints->GetPoint(pts[0], lastPoint);
74

Ken Martin's avatar
Ken Martin committed
75
    for (vtkIdType p = 1; p < npts; ++p)
76
    {
Ken Martin's avatar
Ken Martin committed
77 78 79 80 81 82
      // Create the new cell
      vtkIdType cellId = newLines->InsertNextCell(this->NumberOfSubdivisions);
      output->GetCellData()->CopyData(input->GetCellData(), i, cellId);

      double curPoint[3];
      newPoints->GetPoint(pts[p], curPoint);
83

Ken Martin's avatar
Ken Martin committed
84
      // Find w, a unit vector pointing from the center of the
luz.paz's avatar
luz.paz committed
85
      // earth directly in between the two endpoints.
Ken Martin's avatar
Ken Martin committed
86 87
      double w[3];
      for (int c = 0; c < 3; ++c)
88
      {
Ken Martin's avatar
Ken Martin committed
89
        w[c] = (lastPoint[c] + curPoint[c])/2.0;
90
      }
Ken Martin's avatar
Ken Martin committed
91
      vtkMath::Normalize(w);
92

Ken Martin's avatar
Ken Martin committed
93 94 95 96
      // The center of the circle used to draw the arc is a
      // point along the vector w scaled by the explode factor.
      double center[3];
      for (int c = 0; c < 3; ++c)
97
      {
Ken Martin's avatar
Ken Martin committed
98
        center[c] = this->ExplodeFactor * this->GlobeRadius * w[c];
99
      }
100

Ken Martin's avatar
Ken Martin committed
101 102 103 104 105
      // The vectors u and x are unit vectors pointing from the
      // center of the circle to the two endpoints of the arc,
      // lastPoint and curPoint, respectively.
      double u[3], x[3];
      for (int c = 0; c < 3; ++c)
106
      {
Ken Martin's avatar
Ken Martin committed
107 108
        u[c] = lastPoint[c] - center[c];
        x[c] = curPoint[c] - center[c];
109
      }
Ken Martin's avatar
Ken Martin committed
110 111 112
      double radius = vtkMath::Norm(u);
      vtkMath::Normalize(u);
      vtkMath::Normalize(x);
113

Ken Martin's avatar
Ken Martin committed
114 115
      // Find the angle that the arc spans.
      double theta = acos(vtkMath::Dot(u, x));
116

Ken Martin's avatar
Ken Martin committed
117 118 119 120 121
      // If the vectors u, x point toward the center of the earth, take
      // the larger angle between the vectors.
      // We determine whether u points toward the center of the earth
      // by checking whether the dot product of u and w is negative.
      if (vtkMath::Dot(w, u) < 0)
122
      {
Ken Martin's avatar
Ken Martin committed
123
        theta = 2.0*vtkMath::Pi() - theta;
124
      }
Ken Martin's avatar
Ken Martin committed
125 126 127 128 129 130 131 132 133 134 135 136

      // We need two perpendicular vectors on the plane of the circle
      // in order to draw the circle.  First we calculate n, a vector
      // normal to the circle, by crossing u and w.  Next, we cross
      // n and u in order to get a vector v in the plane of the circle
      // that is perpendicular to u.
      double n[3];
      vtkMath::Cross(u, w, n);
      vtkMath::Normalize(n);
      double v[3];
      vtkMath::Cross(n, u, v);
      vtkMath::Normalize(v);
137

Ken Martin's avatar
Ken Martin committed
138 139 140
      // Use the general equation for a circle in three dimensions
      // to draw an arc from the last point to the current point.
      for (int s = 0; s < this->NumberOfSubdivisions; ++s)
141
      {
Ken Martin's avatar
Ken Martin committed
142 143 144
        double angle = s * theta / (this->NumberOfSubdivisions - 1.0);
        double circlePt[3];
        for (int c = 0; c < 3; ++c)
145
        {
Ken Martin's avatar
Ken Martin committed
146
          circlePt[c] = center[c] + radius*cos(angle)*u[c] + radius*sin(angle)*v[c];
147
        }
Ken Martin's avatar
Ken Martin committed
148 149
        vtkIdType newPt = newPoints->InsertNextPoint(circlePt);
        newLines->InsertCellPoint(newPt);
150
      }
151

Ken Martin's avatar
Ken Martin committed
152
      for (int c = 0; c < 3; ++c)
153
      {
Ken Martin's avatar
Ken Martin committed
154 155 156
        lastPoint[c] = curPoint[c];
      }
    }
157
  }
Ken Martin's avatar
Ken Martin committed
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176

  // Send the data to output.
  output->SetLines(newLines);
  output->SetPoints(newPoints);

  // Clean up.
  newLines->Delete();
  newPoints->Delete();

  return 1;
}

void vtkGeoArcs::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os,indent);
  os << indent << "GlobeRadius: " << this->GlobeRadius << endl;
  os << indent << "ExplodeFactor: " << this->ExplodeFactor << endl;
  os << indent << "NumberOfSubdivisions: " << this->NumberOfSubdivisions << endl;
}