vtkProjectSphereFilter.cxx 19.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
/*=========================================================================

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

=========================================================================*/
#include "vtkProjectSphereFilter.h"

#include "vtkCell.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDoubleArray.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMath.h"
#include "vtkMergePoints.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"

#include <map>

namespace
{
39
  void ConvertXYZToLatLonDepth(double xyz[3], double lonLatDepth[3], double center[3])
40 41
  {
    double dist2 = vtkMath::Distance2BetweenPoints(xyz, center);
42 43 44 45
    lonLatDepth[2] = sqrt(dist2);
    double radianAngle = atan2(xyz[1]-center[1], xyz[0]-center[0]);
    lonLatDepth[0] = radianAngle*180./vtkMath::Pi()-180.;
    lonLatDepth[1] = 90.-acos((xyz[2]-center[2])/lonLatDepth[2])*180./vtkMath::Pi();
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
  }

  template<class data_type>
  void TransformVector(
    double* transformMatrix, data_type* data)
  {
    double d0 = static_cast<double>(data[0]);
    double d1 = static_cast<double>(data[1]);
    double d2 = static_cast<double>(data[2]);
    data[0] = static_cast<data_type>(
      transformMatrix[0]*d0+transformMatrix[1]*d1+transformMatrix[2]*d2);
    data[1] = static_cast<data_type>(
      transformMatrix[3]*d0+transformMatrix[4]*d1+transformMatrix[5]*d2);
    data[2] = static_cast<data_type>(
      transformMatrix[6]*d0+transformMatrix[7]*d1+transformMatrix[8]*d2);
  }
} // end anonymous namespace

vtkStandardNewMacro(vtkProjectSphereFilter);

//-----------------------------------------------------------------------------
vtkProjectSphereFilter::vtkProjectSphereFilter() : SplitLongitude(-180)
{
  this->Center[0] = this->Center[1] = this->Center[2] = 0;
  this->KeepPolePoints = false;
  this->TranslateZ = false;
}

//-----------------------------------------------------------------------------
75
vtkProjectSphereFilter::~vtkProjectSphereFilter() = default;
76 77 78 79 80 81

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

82 83 84
  double center[3];
  this->GetCenter(center);

85
  os << indent << "Center: ("
86 87 88 89 90
     << center[0] << ", "
     << center[1] << ", "
     << center[2] << ")\n";
  os << indent << "KeepPolePoints " << this->GetKeepPolePoints() << "\n";
  os << indent << "TranslateZ " << this->GetTranslateZ() << "\n";
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
}

//-----------------------------------------------------------------------------
int vtkProjectSphereFilter::FillInputPortInformation(int vtkNotUsed(port),
                                                     vtkInformation *info)
{
  info->Remove(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE());
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
  info->Append(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
  return 1;
}

//-----------------------------------------------------------------------------
int vtkProjectSphereFilter::RequestData(vtkInformation *vtkNotUsed(request),
                                        vtkInformationVector **inputVector,
                                        vtkInformationVector *outputVector)
{
  vtkDebugMacro("RequestData");

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation *outInfo = outputVector->GetInformationObject(0);

  vtkPointSet *input
    = vtkPointSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  if(vtkPolyData* poly = vtkPolyData::SafeDownCast(input))
116
  {
117 118 119
    if( poly->GetVerts()->GetNumberOfCells() > 0 ||
        poly->GetLines()->GetNumberOfCells() > 0 ||
        poly->GetStrips()->GetNumberOfCells() > 0 )
120
    {
121 122 123
      vtkErrorMacro("Can only deal with vtkPolyData polys.");
      return 0;
    }
124
  }
125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143

  vtkPointSet *output
    = vtkPointSet::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkNew<vtkIdList> polePointIds;
  this->TransformPointInformation(input, output, polePointIds.GetPointer());
  this->TransformCellInformation(input, output, polePointIds.GetPointer());
  output->GetFieldData()->ShallowCopy(input->GetFieldData());

  vtkDebugMacro("Leaving RequestData");

  return 1;
}

//-----------------------------------------------------------------------------
void vtkProjectSphereFilter::TransformPointInformation(
  vtkPointSet* input, vtkPointSet* output, vtkIdList* polePointIds)
{
  polePointIds->Reset();
144 145 146 147 148

  // Deep copy point data since TransformPointInformation modifies
  // the point data
  output->GetPointData()->DeepCopy(input->GetPointData());

149 150 151 152 153 154 155 156 157
  vtkNew<vtkPoints> points;
  points->SetDataTypeToDouble();
  points->SetNumberOfPoints(input->GetNumberOfPoints());

  double zTranslation = ( this->TranslateZ == true ?
                          this->GetZTranslation(input) : 0. );

  output->SetPoints(points.GetPointer());
  vtkIdType numberOfPoints = input->GetNumberOfPoints();
158
  double minDist2ToCenterLine = VTK_DOUBLE_MAX;
159
  for(vtkIdType i=0;i<numberOfPoints;i++)
160
  {
161 162 163
    double coordIn[3], coordOut[3];
    input->GetPoint(i, coordIn);
    ConvertXYZToLatLonDepth(coordIn, coordOut, this->Center);
164 165 166 167
    // if we allow the user to specify SplitLongitude we have to make
    // sure that we respect their choice since the output of atan
    // is from -180 to 180.
    if(coordOut[0] < this->SplitLongitude)
168
    {
169
      coordOut[0] += 360.;
170
    }
171 172 173 174 175 176 177 178 179 180
    coordOut[2] -= zTranslation;

    // acbauer -- a hack to make the grid look better by forcing it to be flat.
    // leaving this in for now even though it's commented out. if I figure out
    // a proper way to do this i'll replace it.
    // if(this->TranslateZ)
    //   {
    //   coordOut[2] = 0;
    //   }
    points->SetPoint(i, coordOut);
181 182 183 184

    // keep track of the ids of the points that are closest to the
    // centerline between -90 and 90 latitude. this is done as a single
    // pass algorithm.
185 186 187
    double dist2 =
      (coordIn[0]-this->Center[0])*(coordIn[0]-this->Center[0])+
      (coordIn[1]-this->Center[1])*(coordIn[1]-this->Center[1]);
188
    if(dist2 < minDist2ToCenterLine)
189
    {
190 191 192
      // we found a closer point so throw out the previous closest
      // point ids.
      minDist2ToCenterLine = dist2;
193 194
      polePointIds->SetNumberOfIds(1);
      polePointIds->SetId(0, i);
195
    }
196
    else if(dist2 == minDist2ToCenterLine)
197
    {
198 199
      // this point is just as close as the current closest point
      // so we just add it to our list.
200 201
      polePointIds->InsertNextId(i);
    }
202 203
    this->TransformTensors(i, coordIn, output->GetPointData());
  }
204
  this->ComputePointsClosestToCenterLine(minDist2ToCenterLine, polePointIds);
205 206 207 208 209 210 211 212 213 214 215 216 217 218 219
}

//-----------------------------------------------------------------------------
void vtkProjectSphereFilter::TransformCellInformation(
  vtkPointSet* input, vtkPointSet* output, vtkIdList* polePointIds)
{
  // a map from the old point to the newly created point for split cells
  std::map<vtkIdType,vtkIdType> boundaryMap;

  double TOLERANCE = .0001;
  vtkNew<vtkMergePoints> locator;
  locator->InitPointInsertion(output->GetPoints(), output->GetBounds(),
                              output->GetNumberOfPoints());
  double coord[3];
  for(vtkIdType i=0;i<output->GetNumberOfPoints();i++)
220
  {
221 222 223 224 225
    // this is a bit annoying but required for building up the locator properly
    // otherwise it won't either know these points exist or will start
    // counting new points at index 0.
    output->GetPoint(i, coord);
    locator->InsertNextPoint(coord);
226
  }
227 228

  vtkIdType numberOfCells = input->GetNumberOfCells();
229
  vtkCellArray* connectivity = nullptr;
230 231 232
  vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::SafeDownCast(output);
  vtkPolyData* poly = vtkPolyData::SafeDownCast(output);
  if(ugrid)
233
  {
234 235
    ugrid->Allocate(numberOfCells);
    connectivity = ugrid->GetCells();
236
  }
237
  else if(poly)
238
  {
239 240
    poly->Allocate(numberOfCells);
    connectivity = poly->GetPolys();
241
  }
242 243 244 245 246 247 248 249 250 251 252
  output->GetCellData()->CopyAllOn();
  output->GetCellData()->CopyAllocate(input->GetCellData(),
                                      input->GetNumberOfCells());
  vtkPointData* pointData = output->GetPointData();
  pointData->CopyAllOn();
  pointData->CopyAllocate(pointData, output->GetNumberOfPoints());

  vtkNew<vtkIdList> cellPoints;
  vtkNew<vtkIdList> skippedCells;
  vtkIdType mostPointsInCell = 0;
  for(vtkIdType cellId=0;cellId<numberOfCells;cellId++)
253
  {
254 255 256 257 258 259 260 261 262 263 264 265 266
    bool onLeftBoundary = false;
    bool onRightBoundary = false;
    bool leftSideInterior = false; //between SplitLongitude and SplitLongitude+90
    bool rightSideInterior = false;// between SplitLongitude+270 and SplitLongitude+360
    bool middleInterior = false; //between SplitLongitude+90 and SplitLongitude+270

    bool skipCell = false;
    bool splitCell = false;
    double xyz[3];
    input->GetCellPoints(cellId, cellPoints.GetPointer());
    mostPointsInCell = (mostPointsInCell > cellPoints->GetNumberOfIds() ?
                        mostPointsInCell : cellPoints->GetNumberOfIds() );
    for(vtkIdType pt=0;pt<cellPoints->GetNumberOfIds();pt++)
267
    {
268 269
      output->GetPoint(cellPoints->GetId(pt), xyz);
      if(xyz[0] < this->SplitLongitude+TOLERANCE)
270
      {
271
        onLeftBoundary = true;
272
      }
273
      else if(xyz[0] > this->SplitLongitude+360.-TOLERANCE)
274
      {
275
        onRightBoundary = true;
276
      }
277
      else if(xyz[0] < this->SplitLongitude+90.)
278
      {
279
        leftSideInterior = true;
280
      }
281
      else if(xyz[0] > this->SplitLongitude+270.)
282
      {
283
        rightSideInterior = true;
284
      }
285
      else
286
      {
287
        middleInterior = true;
288
      }
289
      if(polePointIds->IsId(cellPoints->GetId(pt)) != -1 && this->KeepPolePoints == false)
290
      {
291 292 293 294
        skipCell = true;
        skippedCells->InsertNextId(cellId);
        continue;
      }
295
    }
296
    if(skipCell)
297
    {
298
      continue;
299
    }
300
    if( (onLeftBoundary || onRightBoundary ) && rightSideInterior && leftSideInterior)
301
    { // this cell stretches across the split longitude
302
      splitCell = true;
303
    }
304
    else if( onLeftBoundary && rightSideInterior )
305
    {
306
      for(vtkIdType pt=0;pt<cellPoints->GetNumberOfIds();pt++)
307
      {
308 309
        output->GetPoint(cellPoints->GetId(pt), xyz);
        if(xyz[0] < this->SplitLongitude+TOLERANCE)
310
        {
311 312 313
          std::map<vtkIdType,vtkIdType>::iterator it =
            boundaryMap.find(cellPoints->GetId(pt));
          if(it == boundaryMap.end())
314
          { // need to create another point
315 316 317 318 319
            xyz[0] += 360.;
            vtkIdType id = locator->InsertNextPoint(xyz);
            boundaryMap[cellPoints->GetId(pt)] = id;
            pointData->CopyData(pointData, cellPoints->GetId(pt), id);
            cellPoints->SetId(pt, id);
320
          }
321
          else
322
          {
323 324 325 326
            cellPoints->SetId(pt, it->second);
          }
        }
      }
327
    }
328
    else if( onRightBoundary && leftSideInterior )
329
    {
330
      for(vtkIdType pt=0;pt<cellPoints->GetNumberOfIds();pt++)
331
      {
332 333
        output->GetPoint(cellPoints->GetId(pt), xyz);
        if(xyz[0] > this->SplitLongitude+360.-TOLERANCE)
334
        {
335 336 337
          std::map<vtkIdType,vtkIdType>::iterator it =
            boundaryMap.find(cellPoints->GetId(pt));
          if(it == boundaryMap.end())
338
          { // need to create another point
339 340 341 342 343
            xyz[0] -= 360.;
            vtkIdType id = locator->InsertNextPoint(xyz);
            boundaryMap[cellPoints->GetId(pt)] = id;
            pointData->CopyData(pointData, cellPoints->GetId(pt), id);
            cellPoints->SetId(pt, id);
344
          }
345
          else
346
          {
347 348 349 350
            cellPoints->SetId(pt, it->second);
          }
        }
      }
351
    }
352
    else if( (onLeftBoundary || onRightBoundary ) && middleInterior)
353
    {
354
      splitCell = true;
355
    }
356
    else if( leftSideInterior && rightSideInterior)
357
    {
358
      splitCell = true;
359
    }
360
    if(splitCell)
361
    {
362 363
      this->SplitCell(input, output, cellId, locator.GetPointer(), connectivity, 0);
      this->SplitCell(input, output, cellId, locator.GetPointer(), connectivity, 1);
364
    }
365
    else if(ugrid)
366
    {
367 368 369
      ugrid->InsertNextCell(input->GetCellType(cellId), cellPoints.GetPointer());
      output->GetCellData()->CopyData(input->GetCellData(), cellId,
                                      output->GetNumberOfCells()-1);
370
    }
371
    else if(poly)
372
    {
373 374 375 376
      poly->InsertNextCell(input->GetCellType(cellId), cellPoints.GetPointer());
      output->GetCellData()->CopyData(input->GetCellData(), cellId,
                                      output->GetNumberOfCells()-1);
    }
377
  }
378 379

  if(poly)
380
  {
381 382 383 384
    // we have to rebuild the polydata cell data structures since when
    // we split a cell we don't do it right away due to the expense
    poly->DeleteCells();
    poly->BuildCells();
385
  }
386 387 388 389 390

  // deal with cell data
  std::vector<double> weights(mostPointsInCell);
  vtkIdType skipCounter = 0;
  for(vtkIdType cellId=0;cellId<input->GetNumberOfCells();cellId++)
391
  {
392
    if(skippedCells->IsId(cellId) != -1)
393
    {
394 395 396
      skippedCells->DeleteId(cellId);
      skipCounter++;
      continue;
397
    }
398 399 400 401 402 403
    int subId = 0;
    double parametricCenter[3];
    vtkCell* cell = input->GetCell(cellId);
    cell->GetParametricCenter(parametricCenter);
    cell->EvaluateLocation(subId, parametricCenter, coord, &weights[0]);
    this->TransformTensors(cellId-skipCounter, coord, output->GetCellData());
404
  }
405 406 407 408 409
}

//-----------------------------------------------------------------------------
void vtkProjectSphereFilter::TransformTensors(
  vtkIdType pointId, double* coord, vtkDataSetAttributes* dataArrays  )
410
{
411 412 413 414 415 416 417 418 419 420 421 422 423 424
   double theta = atan2(
     sqrt( (coord[0]-this->Center[0])*(coord[0]-this->Center[0]) +
           (coord[1]-this->Center[1])*(coord[1]-this->Center[1]) ),
     coord[2]-this->Center[2] );
   double phi = atan2( coord[1]-this->Center[1], coord[0]-this->Center[0] );
   double sinTheta = sin(theta);
   double cosTheta = cos(theta);
   double sinPhi = sin(phi);
   double cosPhi = cos(phi);
   double transformMatrix[9] = {
     -sinPhi, cosPhi, 0.,
     cosTheta*cosPhi, cosTheta*sinPhi, -sinTheta,
     sinTheta*cosPhi, sinTheta*sinPhi, cosTheta };
   for(int i=0;i<dataArrays->GetNumberOfArrays();i++)
425
   {
426 427
     vtkDataArray* array = dataArrays->GetArray(i);
     if(array->GetNumberOfComponents() == 3)
428
     {
429
       switch (array->GetDataType())
430
       {
431 432 433 434 435
         vtkTemplateMacro(TransformVector(
                            transformMatrix,
                            static_cast<VTK_TT *>(array->GetVoidPointer(pointId*array->GetNumberOfComponents())) ) );
       }
     }
436 437
   }
}
438 439 440 441 442 443 444

//-----------------------------------------------------------------------------
double vtkProjectSphereFilter::GetZTranslation(vtkPointSet* input)
{
  double maxRadius2 = 0; // squared radius
  double coord[3];
  for(vtkIdType i=0;i<input->GetNumberOfPoints();i++)
445
  {
446 447 448
    input->GetPoint(i, coord);
    double dist2 = vtkMath::Distance2BetweenPoints(coord, this->Center);
    if(dist2 > maxRadius2)
449
    {
450 451
      maxRadius2 = dist2;
    }
452
  }
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468
  return sqrt(maxRadius2);
}

//-----------------------------------------------------------------------------
void vtkProjectSphereFilter::SplitCell(
  vtkPointSet* input, vtkPointSet* output, vtkIdType inputCellId,
  vtkIncrementalPointLocator* locator, vtkCellArray* connectivity,
  int splitSide)
{
  // i screw up the canonical ordering of the cell but apparently this
  // gets fixed by vtkCell::Clip().
  vtkCell* cell = input->GetCell(inputCellId);
  vtkNew<vtkDoubleArray> cellScalars;
  cellScalars->SetNumberOfTuples(cell->GetNumberOfPoints());
  double coord[3];
  for(vtkIdType pt=0;pt<cell->GetNumberOfPoints();pt++)
469
  {
470 471
    output->GetPoint(cell->GetPointId(pt), coord);
    if(splitSide == 0 && coord[0] > this->SplitLongitude+180.)
472
    {
473
      coord[0] -= 360.;
474
    }
475
    else if(splitSide == 1 && coord[0] < this->SplitLongitude+180.)
476
    {
477
      coord[0] += 360.;
478
    }
479 480
    cellScalars->SetValue(pt, coord[0]);
    cell->GetPoints()->SetPoint(pt, coord);
481
  }
482 483 484 485 486 487 488 489
  vtkIdType numberOfCells = output->GetNumberOfCells();
  double splitLocation = (splitSide == 0 ? -180. : 180.);
  cell->Clip(splitLocation, cellScalars.GetPointer(), locator, connectivity,
             output->GetPointData(), output->GetPointData(), input->GetCellData(),
             inputCellId, output->GetCellData(), splitSide);
  // if the grid was an unstructured grid we have to update the cell
  // types and locations for the created cells.
  if(vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::SafeDownCast(output))
490
  {
491
    this->SetCellInformation(ugrid, cell, output->GetNumberOfCells()-numberOfCells);
492
  }
493 494 495 496 497 498 499
}

//-----------------------------------------------------------------------------
void vtkProjectSphereFilter::SetCellInformation(
  vtkUnstructuredGrid* output, vtkCell* cell, vtkIdType numberOfNewCells)
{
  for(vtkIdType i=0;i<numberOfNewCells;i++)
500
  {
501 502 503 504 505 506 507 508 509
    vtkIdType prevCellId = output->GetNumberOfCells()+i-numberOfNewCells-1;
    vtkIdType newCellId = prevCellId + 1;
    vtkIdType *pts, numPts;
    vtkIdType loc = output->GetCellLocationsArray()->GetValue(prevCellId);
    output->GetCells()->GetCell(loc,numPts,pts);

    output->GetCellLocationsArray()->InsertNextValue(loc+numPts+1);
    output->GetCells()->GetCell(loc+numPts+1,numPts,pts);
    if(cell->GetCellDimension() == 0)
510
    {
511
      if(numPts > 2)
512
      {
513
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_POLY_VERTEX);
514
      }
515
      else
516
      {
517 518
        vtkErrorMacro("Cannot handle 0D cell with " << numPts << " number of points.");
      }
519
    }
520
    else if(cell->GetCellDimension() == 1)
521
    {
522
      if(numPts == 2)
523
      {
524
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_LINE);
525
      }
526
      else if(numPts > 2)
527
      {
528
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_POLY_LINE);
529
      }
530
      else
531
      {
532 533
        vtkErrorMacro("Cannot handle 1D cell with " << numPts << " number of points.");
      }
534
    }
535
    else if(cell->GetCellDimension() == 2)
536
    {
537
      if(numPts == 3)
538
      {
539
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_TRIANGLE);
540
      }
541
      else if(numPts > 3 && cell->GetCellType() == VTK_TRIANGLE_STRIP)
542
      {
543
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_TRIANGLE_STRIP);
544
      }
545
      else if(numPts == 4)
546
      {
547
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_QUAD);
548
      }
549
      else
550
      {
551 552
        vtkErrorMacro("Cannot handle 2D cell with " << numPts << " number of points.");
      }
553
    }
554
    else  // 3D cell
555
    {
556
      if(numPts == 4)
557
      {
558
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_TETRA);
559
      }
560
      else if(numPts == 5)
561
      {
562
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_PYRAMID);
563
      }
564
      else if(numPts == 6)
565
      {
566
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_WEDGE);
567
      }
568
      else if(numPts == 8)
569
      {
570
        output->GetCellTypesArray()->InsertValue(newCellId, VTK_HEXAHEDRON);
571
      }
572
      else
573
      {
574 575 576
        vtkErrorMacro("Unknown 3D cell type.");
      }
    }
577
  }
578
}