TestVTKMDataSet.cxx 12.7 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
#include "vtkmDataSet.h"

#include "vtkCell.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
#include "vtkIdTypeArray.h"
#include "vtkImageData.h"
#include "vtkImageDataToPointSet.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkStructuredGrid.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"
19
#include "vtkmDataArray.h"
20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

#include <vtkm/cont/ArrayHandleUniformPointCoordinates.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>

#include <array>
#include <random>

namespace
{

//-----------------------------------------------------------------------------
class TestError
{
public:
  TestError(const std::string& message, int line)
35 36
    : Message(message)
    , Line(line)
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
  {
  }

  void PrintMessage(std::ostream& out) const
  {
    out << "Error at line " << this->Line << ": " << this->Message << "\n";
  }

private:
  std::string Message;
  int Line;
};

#define RAISE_TEST_ERROR(msg) throw TestError((msg), __LINE__)

52 53 54
#define TEST_VERIFY(cond, msg)                                                                     \
  if (!(cond))                                                                                     \
  RAISE_TEST_ERROR((msg))
55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73

inline bool IsEqualFloat(double a, double b, double e = 1e-6f)
{
  return (std::abs(a - b) <= e);
}

//-----------------------------------------------------------------------------
inline void TestEqualCells(vtkCell* c1, vtkCell* c2)
{
  TEST_VERIFY(c1->GetCellType() == c2->GetCellType(), "Cell types don't match");
  TEST_VERIFY(c1->GetNumberOfPoints() == c2->GetNumberOfPoints(), "Cell sizes don't match");
  for (vtkIdType i = 0; i < c1->GetNumberOfPoints(); ++i)
  {
    TEST_VERIFY(c1->GetPointId(i) == c2->GetPointId(i), "Cell point-ids don't match");
  }
}

inline void TestEqualVtkArrays(vtkAbstractArray* a1, vtkAbstractArray* a2)
{
74 75 76 77
  TEST_VERIFY(std::string(a1->GetName()) == std::string(a2->GetName()), "Array names don't match");
  TEST_VERIFY(a1->GetDataType() == a2->GetDataType(), "Array data-types don't match");
  TEST_VERIFY(
    a1->GetNumberOfTuples() == a2->GetNumberOfTuples(), "Array number of tuples don't match");
78
  TEST_VERIFY(a1->GetNumberOfComponents() == a2->GetNumberOfComponents(),
79
    "Array number of components don't match");
80 81 82 83 84 85 86 87 88 89

  auto da1 = vtkDataArray::SafeDownCast(a1);
  auto da2 = vtkDataArray::SafeDownCast(a2);
  double range1[2], range2[2];
  int numComps = da1->GetNumberOfComponents();
  for (int i = 0; i < numComps; ++i)
  {
    da1->GetRange(range1, i);
    da2->GetRange(range2, i);
    TEST_VERIFY(IsEqualFloat(range1[0], range2[0]) && IsEqualFloat(range1[1], range2[1]),
90
      "Array ranges don't match");
91 92 93 94 95
  }
}

void TestDataSets(vtkDataSet* dsVtk, vtkDataSet* dsVtkm)
{
96 97 98 99
  TEST_VERIFY(
    dsVtk->GetNumberOfPoints() == dsVtkm->GetNumberOfPoints(), "Number of points don't match");
  TEST_VERIFY(
    dsVtk->GetNumberOfCells() == dsVtkm->GetNumberOfCells(), "Number of cells don't match");
100 101 102 103
  double bounds1[6], bounds2[6];
  dsVtk->GetBounds(bounds1);
  dsVtkm->GetBounds(bounds2);
  TEST_VERIFY(IsEqualFloat(bounds1[0], bounds2[0]) && IsEqualFloat(bounds1[1], bounds2[1]) &&
104 105 106
      IsEqualFloat(bounds1[2], bounds2[2]) && IsEqualFloat(bounds1[3], bounds2[3]) &&
      IsEqualFloat(bounds1[4], bounds2[4]) && IsEqualFloat(bounds1[5], bounds2[5]),
    "Bounds don't match");
107 108 109 110 111 112 113

  vtkIdType numPoints = dsVtk->GetNumberOfPoints();
  for (vtkIdType i = 0; i < numPoints; ++i)
  {
    double x1[3], x2[3];
    dsVtk->GetPoint(i, x1);
    dsVtkm->GetPoint(i, x2);
114 115 116
    TEST_VERIFY(
      IsEqualFloat(x1[0], x2[0]) && IsEqualFloat(x1[1], x2[1]) && IsEqualFloat(x1[2], x2[2]),
      "'GetPoint` results don't match");
117 118 119 120 121 122 123

    vtkNew<vtkIdList> cellIds1, cellIds2;
    dsVtk->GetPointCells(i, cellIds1);
    dsVtkm->GetPointCells(i, cellIds2);
    cellIds1->Sort();
    cellIds2->Sort();
    TEST_VERIFY(cellIds1->GetNumberOfIds() == cellIds2->GetNumberOfIds(),
124
      "`GetPointCells` results don't match");
125 126
    for (vtkIdType j = 0; j < cellIds1->GetNumberOfIds(); ++j)
    {
127
      TEST_VERIFY(cellIds1->GetId(j) == cellIds2->GetId(j), "`GetPointCells` results don't match");
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144
    }
  }

  vtkIdType numCells = dsVtk->GetNumberOfCells();
  for (vtkIdType i = 0; i < numCells; ++i)
  {
    TestEqualCells(dsVtk->GetCell(i), dsVtkm->GetCell(i));

    vtkNew<vtkGenericCell> gc1, gc2;
    dsVtk->GetCell(i, gc1);
    dsVtkm->GetCell(i, gc2);
    TestEqualCells(gc1, gc2);

    double bds1[6], bds2[6];
    dsVtk->GetCellBounds(i, bds1);
    dsVtkm->GetCellBounds(i, bds2);
    TEST_VERIFY(IsEqualFloat(bds1[0], bds2[0]) && IsEqualFloat(bds1[1], bds2[1]) &&
145 146 147
        IsEqualFloat(bds1[2], bds2[2]) && IsEqualFloat(bds1[3], bds2[3]) &&
        IsEqualFloat(bds1[4], bds2[4]) && IsEqualFloat(bds1[5], bds2[5]),
      "Cell bounds don't match");
148

149
    TEST_VERIFY(dsVtk->GetCellType(i) == dsVtkm->GetCellType(i), "Cell types don't match");
150 151 152 153 154 155

    vtkNew<vtkIdList> ptIds1, ptIds2;
    dsVtk->GetCellPoints(i, ptIds1);
    dsVtkm->GetCellPoints(i, ptIds2);
    for (vtkIdType j = 0; j < ptIds1->GetNumberOfIds(); ++j)
    {
156
      TEST_VERIFY(ptIds1->GetId(j) == ptIds2->GetId(j), "`GetCellPoints` results don't match");
157 158 159 160
    }
  }

  std::default_random_engine engine;
161 162
  std::uniform_real_distribution<double> d1(bounds1[0], bounds1[1]), d2(bounds1[2], bounds1[3]),
    d3(bounds1[4], bounds1[5]);
163 164 165 166 167 168 169 170 171 172 173 174 175
  static constexpr int numSamples = 100;
  for (int i = 0; i < numSamples; ++i)
  {
    double x[3] = { d1(engine), d2(engine), d3(engine) };

    auto pid1 = dsVtk->FindPoint(x);
    auto pid2 = dsVtkm->FindPoint(x);
    if (pid1 != pid2)
    {
      TEST_VERIFY(pid1 != -1 && pid2 != -1, "`FindPoint` results don't match");
      double x1[3], x2[3];
      dsVtk->GetPoint(pid1, x1);
      dsVtkm->GetPoint(pid1, x2);
176 177 178
      TEST_VERIFY(IsEqualFloat(
                    vtkMath::Distance2BetweenPoints(x, x1), vtkMath::Distance2BetweenPoints(x, x2)),
        "`FindPoint` results don't match");
179 180 181 182 183 184 185 186
    }

    int subId = 0;
    double pcoords1[3], pcoords2[3];
    double weights1[8], weights2[8];
    auto cid1 = dsVtk->FindCell(x, nullptr, -1, 1e-6, subId, pcoords1, weights1);
    auto cid2 = dsVtkm->FindCell(x, nullptr, -1, 1e-6, subId, pcoords2, weights2);

187 188 189 190
    // vtkDataSet and vtkmDataSet may find different cells if the point is too close to
    // the boundary of those cells
    if (cid1 != cid2)
    {
191
      if (cid2 >= 0)
192
      {
193 194 195 196 197 198
        // check if the point is inside or close to the vtkmDataSet found cell
        vtkCell* cell = dsVtk->GetCell(cid2);
        double dist2 = 0, pcoords[3], weights[8];
        if (cell->EvaluatePosition(x, nullptr, subId, pcoords, dist2, weights) == 0) // outside?
        {
          TEST_VERIFY(IsEqualFloat(cell->GetParametricDistance(pcoords), 0.0, 1e-3),
199
            "`FindCell` incorrect result by vtkmDataSet");
200
        }
201 202 203
      }
    }
    else if (cid1 == -1)
204 205 206
    {
      continue;
    }
207
    else
208
    {
209
      TEST_VERIFY(IsEqualFloat(pcoords1[0], pcoords2[0]) &&
210 211
          IsEqualFloat(pcoords1[1], pcoords2[1]) && IsEqualFloat(pcoords1[2], pcoords2[2]),
        "`FindCell` pcoords don't match");
212 213 214
      int count = dsVtk->GetCell(cid1)->GetNumberOfPoints();
      for (int j = 0; j < count; ++j)
      {
215
        TEST_VERIFY(IsEqualFloat(weights1[j], weights2[j]), "`FindCell` weights don't match");
216
      }
217 218 219 220 221 222
    }
  }

  // test fields
  int numPointFields = dsVtk->GetPointData()->GetNumberOfArrays();
  TEST_VERIFY(numPointFields == dsVtkm->GetPointData()->GetNumberOfArrays(),
223
    "Number of point-fields don't match");
224 225
  for (int i = 0; i < numPointFields; ++i)
  {
226
    TestEqualVtkArrays(dsVtk->GetPointData()->GetArray(i), dsVtkm->GetPointData()->GetArray(i));
227 228 229 230
  }

  int numCellFields = dsVtk->GetCellData()->GetNumberOfArrays();
  TEST_VERIFY(numCellFields == dsVtkm->GetCellData()->GetNumberOfArrays(),
231
    "Number of cell-fields don't match");
232 233
  for (int i = 0; i < numCellFields; ++i)
  {
234
    TestEqualVtkArrays(dsVtk->GetCellData()->GetArray(i), dsVtkm->GetCellData()->GetArray(i));
235 236 237 238
  }
}

//-----------------------------------------------------------------------------
239
inline void CoordsCopy(const vtkm::cont::CoordinateSystem& coords, vtkPoints* points)
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
{
  auto ptsPortal = coords.GetData().GetPortalConstControl();
  auto numPoints = coords.GetNumberOfPoints();

  points->SetDataTypeToFloat();
  points->SetNumberOfPoints(numPoints);
  auto ptsArray = vtkFloatArray::SafeDownCast(points->GetData());
  for (vtkIdType i = 0; i < numPoints; ++i)
  {
    auto pt = ptsPortal.Get(i);
    float tuple[3] = { pt[0], pt[1], pt[2] };
    ptsArray->SetTypedTuple(i, tuple);
  }
}

255 256
inline void FieldCopy(
  const vtkm::cont::ArrayHandle<float>& src, const char* name, vtkFloatArray* dst)
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283
{
  auto portal = src.GetPortalConstControl();
  vtkm::Id length = portal.GetNumberOfValues();

  dst->SetName(name);
  dst->SetNumberOfComponents(1);
  dst->SetNumberOfTuples(length);
  for (vtkm::Id i = 0; i < length; ++i)
  {
    dst->SetValue(i, portal.Get(i));
  }
}

//-----------------------------------------------------------------------------
vtkm::cont::testing::MakeTestDataSet Maker;

void TestUniformDataSet()
{
  auto dataset = Maker.Make3DUniformDataSet0();
  auto coords =
    dataset.GetCoordinateSystem().GetData().Cast<vtkm::cont::ArrayHandleUniformPointCoordinates>();
  auto portal = coords.GetPortalConstControl();
  auto dims = portal.GetDimensions();
  auto origin = portal.GetOrigin();
  auto spacing = portal.GetSpacing();

  vtkNew<vtkFloatArray> pointField, cellField;
284 285 286 287
  FieldCopy(dataset.GetField("pointvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "pointvar", pointField);
  FieldCopy(dataset.GetField("cellvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "cellvar", cellField);
288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310

  vtkNew<vtkImageData> imageData;
  imageData->SetDimensions(dims[0], dims[1], dims[2]);
  imageData->SetOrigin(origin[0], origin[1], origin[2]);
  imageData->SetSpacing(spacing[0], spacing[1], spacing[2]);
  imageData->GetPointData()->AddArray(pointField);
  imageData->GetCellData()->AddArray(cellField);

  vtkNew<vtkImageDataToPointSet> voxToHex;
  voxToHex->SetInputData(imageData);
  voxToHex->Update();

  auto dsVtk = voxToHex->GetOutput();

  vtkNew<vtkmDataSet> dsVtkm;
  dsVtkm->SetVtkmDataSet(dataset);

  TestDataSets(dsVtk, dsVtkm);
}

void TestCurvilinearDataSet()
{
  auto dataset = Maker.Make3DRegularDataSet0();
311
  auto dims = dataset.GetCellSet().Cast<vtkm::cont::CellSetStructured<3> >().GetPointDimensions();
312 313 314 315 316

  vtkNew<vtkPoints> points;
  CoordsCopy(dataset.GetCoordinateSystem(), points);

  vtkNew<vtkFloatArray> pointField, cellField;
317 318 319 320
  FieldCopy(dataset.GetField("pointvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "pointvar", pointField);
  FieldCopy(dataset.GetField("cellvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "cellvar", cellField);
321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358

  vtkNew<vtkStructuredGrid> dsVtk;
  dsVtk->SetDimensions(dims[0], dims[1], dims[2]);
  dsVtk->SetPoints(points);
  dsVtk->GetPointData()->AddArray(pointField);
  dsVtk->GetCellData()->AddArray(cellField);

  vtkNew<vtkmDataSet> dsVtkm;
  dsVtkm->SetVtkmDataSet(dataset);

  TestDataSets(dsVtk, dsVtkm);
}

void TestExplicitDataSet()
{
  auto dataset = Maker.Make3DExplicitDataSetZoo();

  vtkNew<vtkPoints> points;
  CoordsCopy(dataset.GetCoordinateSystem(), points);

  auto& cellset = *dataset.GetCellSet().GetCellSetBase();
  vtkm::Id numCells = cellset.GetNumberOfCells();

  vtkNew<vtkUnsignedCharArray> shapes;
  vtkNew<vtkCellArray> connectivity;
  shapes->SetNumberOfComponents(1);
  shapes->SetNumberOfTuples(numCells);
  for (vtkm::Id i = 0; i < numCells; ++i)
  {
    shapes->SetValue(i, cellset.GetCellShape(i));

    vtkIdType ptIds[8];
    int count = cellset.GetNumberOfPointsInCell(i);
    cellset.GetCellPointIds(i, ptIds);
    connectivity->InsertNextCell(count, ptIds);
  }

  vtkNew<vtkFloatArray> pointField, cellField;
359 360 361 362
  FieldCopy(dataset.GetField("pointvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "pointvar", pointField);
  FieldCopy(dataset.GetField("cellvar").GetData().Cast<vtkm::cont::ArrayHandle<float> >(),
    "cellvar", cellField);
363 364 365

  vtkNew<vtkUnstructuredGrid> dsVtk;
  dsVtk->SetPoints(points);
366
  dsVtk->SetCells(shapes, connectivity);
367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400
  dsVtk->GetPointData()->AddArray(pointField);
  dsVtk->GetCellData()->AddArray(cellField);

  vtkNew<vtkmDataSet> dsVtkm;
  dsVtkm->SetVtkmDataSet(dataset);

  TestDataSets(dsVtk, dsVtkm);
}

} // anonymous namespace

//-----------------------------------------------------------------------------
int TestVTKMDataSet(int, char*[])
try
{
  std::cout << "Testing Uniform DataSet\n";
  TestUniformDataSet();
  std::cout << "Passed\n";

  std::cout << "Testing Curvilinear DataSet\n";
  TestCurvilinearDataSet();
  std::cout << "Passed\n";

  std::cout << "Testing Explicit DataSet\n";
  TestExplicitDataSet();
  std::cout << "Passed\n";

  return EXIT_SUCCESS;
}
catch (const TestError& e)
{
  e.PrintMessage(std::cout);
  return 1;
}