Testing/Cxx/BoxClipTriangulate.cxx

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

  Program:   Visualization Toolkit
  Module:    BoxClipTriangulate.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 (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

// This test makes sure that vtkBoxClipDataSet correctly triangulates all cell
// types.

#include "vtkBoxClipDataSet.h"
#include "vtkCellArray.h"
#include "vtkDataSetSurfaceFilter.h"
#include "vtkDoubleArray.h"
#include "vtkMath.h"
#include "vtkPoints.h"
#include "vtkPolyData.h"
#include "vtkTriangle.h"
#include "vtkUnstructuredGrid.h"

#include "vtkExtractEdges.h"
#include "vtkPolyDataMapper.h"
#include "vtkActor.h"
#include "vtkRenderer.h"
#include "vtkRenderWindow.h"
#include "vtkRenderWindowInteractor.h"

#include "vtkSmartPointer.h"
#define VTK_CREATE(type, var) \
  vtkSmartPointer<type> var = vtkSmartPointer<type>::New()

#include <time.h>

#include <vector>

const int NumPoints = 13;
const double PointData[NumPoints*3] = {
  0.0, 0.0, 0.0,
  0.0, 1.0, 0.0,
  1.0, 0.0, 0.0,
  1.0, 1.0, 0.0,
  2.0, 0.0, 0.0,
  2.0, 1.0, 0.0,

  0.0, 0.0, 1.0,
  0.0, 1.0, 1.0,
  1.0, 0.0, 1.0,
  1.0, 1.0, 1.0,
  2.0, 0.0, 1.0,
  2.0, 1.0, 1.0,
  2.0, 0.5, 1.0
};


const vtkIdType NumTriStripCells = 1;
const vtkIdType TriStripCells[] = {
  6, 1, 0, 3, 2, 5, 4
};

const vtkIdType NumQuadCells = 2;
const vtkIdType QuadCells[] = {
  4, 0, 2, 3, 1,
  4, 2, 4, 5, 3
};

const vtkIdType NumPixelCells = 2;
const vtkIdType PixelCells[] = {
  4, 0, 2, 1, 3,
  4, 2, 4, 3, 5
};

const vtkIdType NumPolyCells = 3;
const vtkIdType PolyCells[] = {
  4, 0, 2, 3, 1,
  3, 2, 4, 5,
  5, 6, 8, 12, 9, 7
};


const vtkIdType NumHexCells = 2;
const vtkIdType HexCells[] = {
  8, 6, 8, 2, 0, 7, 9, 3, 1,
  8, 4, 2, 8, 10, 5, 3, 9, 11
};
const vtkIdType NumExpectedHexSurfacePolys = 20;

const vtkIdType NumVoxelCells = 2;
const vtkIdType VoxelCells[] = {
  8, 0, 2, 1, 3, 6, 8, 7, 9,
  8, 10, 8, 11, 9, 4, 2, 5, 3
};
const vtkIdType NumExpectedVoxelSurfacePolys = 20;

const vtkIdType NumWedgeCells = 4;
const vtkIdType WedgeCells[] = {
  6, 0, 1, 2, 6, 7, 8,
  6, 7, 8, 9, 1, 2, 3,
  6, 8, 11, 9, 2, 5, 3,
  6, 2, 5, 4, 8, 11, 10
};
const vtkIdType NumExpectedWedgeSurfacePolys = 20;

const vtkIdType NumPyramidCells = 2;
const vtkIdType PyramidCells[] = {
  5, 8, 9, 3, 2, 0,
  5, 2, 3, 9, 8, 12
};
const vtkIdType NumExpectedPyramidSurfacePolys = 8;

class BoxClipTriangulateFailed { };

//-----------------------------------------------------------------------------

static void CheckWinding(vtkBoxClipDataSet* alg)
{
  alg->Update();
  vtkUnstructuredGrid* data = alg->GetOutput();

  vtkPoints *points = data->GetPoints();

  vtkCellArray *cells = data->GetCells();
  cells->InitTraversal();

  vtkIdType npts, *pts;
  while (cells->GetNextCell(npts, pts))
    {
    if (npts != 4)
      {
      std::cout << "Weird.  I got something that is not a tetrahedra." << std::endl;
      continue;
      }

    double p0[3], p1[3], p2[3], p3[3];
    points->GetPoint(pts[0], p0);
    points->GetPoint(pts[1], p1);
    points->GetPoint(pts[2], p2);
    points->GetPoint(pts[3], p3);

    // If the winding is correct, the normal to triangle p0,p1,p2 should point
    // towards p3.
    double v0[3], v1[3];
    v0[0] = p1[0] - p0[0];    v0[1] = p1[1] - p0[1];    v0[2] = p1[2] - p0[2];
    v1[0] = p2[0] - p0[0];    v1[1] = p2[1] - p0[1];    v1[2] = p2[2] - p0[2];

    double n[3];
    vtkMath::Cross(v0, v1, n);

    double d[3];
    d[0] = p3[0] - p0[0];    d[1] = p3[1] - p0[1];    d[2] = p3[2] - p0[2];

    if (vtkMath::Dot(n, d) < 0)
      {
      std::cout << "Found a tetrahedra with bad winding." << std::endl;
      throw BoxClipTriangulateFailed();
      }
    }
}

//-----------------------------------------------------------------------------

static vtkSmartPointer<vtkUnstructuredGrid> BuildInput(int type,
                                                       vtkIdType numcells,
                                                       const vtkIdType *cells)
{
  vtkIdType i;

  VTK_CREATE(vtkUnstructuredGrid, input);

  // Randomly shuffle the points to possibly test various tessellations.
  // Make a map from original point orderings to new point orderings.
  std::vector<vtkIdType> idMap;
  std::vector<vtkIdType> idsLeft;

  for (i = 0; i < NumPoints; i++)
    {
    idsLeft.push_back(i);
    }

  while (!idsLeft.empty())
    {
    vtkIdType next
      = vtkMath::Round(vtkMath::Random(-0.49, idsLeft.size() - 0.51));
    std::vector<vtkIdType>::iterator nextp = idsLeft.begin() + next;
    idMap.push_back(*nextp);
    idsLeft.erase(nextp, nextp + 1);
    }

  // Build shuffled points.
  VTK_CREATE(vtkPoints, points);
  points->SetNumberOfPoints(NumPoints);
  for (i = 0; i < NumPoints; i++)
    {
    points->SetPoint(idMap[i], PointData + 3*i);
    }
  input->SetPoints(points);

  // Add the cells with indices properly mapped.
  VTK_CREATE(vtkIdList, ptIds);
  const vtkIdType *c = cells;
  for (i = 0; i < numcells; i++)
    {
    vtkIdType npts = *c;  c++;
    ptIds->Initialize();
    for (vtkIdType j = 0; j < npts; j++)
      {
      ptIds->InsertNextId(idMap[*c]);
      c++;
      }
    input->InsertNextCell(type, ptIds);
    }

  return input;
}

//-----------------------------------------------------------------------------

static void Check2DPrimitive(int type, vtkIdType numcells,
                             const vtkIdType *cells)
{
  vtkSmartPointer<vtkUnstructuredGrid> input
    = BuildInput(type, numcells, cells);

  VTK_CREATE(vtkBoxClipDataSet, clipper);
  clipper->SetInputData(input);
  // Clip nothing.
  clipper->SetBoxClip(0.0, 2.0, 0.0, 1.0, 0.0, 1.0);
  clipper->Update();

  vtkUnstructuredGrid *output = clipper->GetOutput();

  if (output->GetNumberOfCells() < 1)
    {
    std::cout << "Output has no cells!" << std::endl;
    throw BoxClipTriangulateFailed();
    }

  // Check to make sure all the normals point in the z direction.
  vtkCellArray *outCells = output->GetCells();
  outCells->InitTraversal();
  vtkIdType npts, *pts;
  while (outCells->GetNextCell(npts, pts))
    {
    if (npts != 3)
      {
      std::cout << "Got a primitive that is not a triangle!" << std::endl;
      throw BoxClipTriangulateFailed();
      }

    double n[3];
    vtkTriangle::ComputeNormal(output->GetPoints(), npts, pts, n);
    if ((n[0] > 0.1) || (n[1] > 0.1) || (n[2] < 0.9))
      {
      std::cout << "Primitive is facing the wrong way!" << std::endl;
      throw BoxClipTriangulateFailed();
      }
    }
}

//-----------------------------------------------------------------------------

static void Check3DPrimitive(int type, vtkIdType numcells,
                             const vtkIdType *cells, vtkIdType numSurfacePolys)
{
  vtkSmartPointer<vtkUnstructuredGrid> input
    = BuildInput(type, numcells, cells);

  VTK_CREATE(vtkBoxClipDataSet, clipper);
  clipper->SetInputData(input);
  // Clip nothing.
  clipper->SetBoxClip(0.0, 2.0, 0.0, 1.0, 0.0, 1.0);
  clipper->Update();

  vtkUnstructuredGrid *output = clipper->GetOutput();

  if (output->GetNumberOfCells() < 1)
    {
    std::cout << "Output has no cells!" << std::endl;
    throw BoxClipTriangulateFailed();
    }

#if 0
  VTK_CREATE(vtkExtractEdges, edges);
  edges->SetInput(output);
  VTK_CREATE(vtkPolyDataMapper, mapper);
  mapper->SetInputConnection(0, edges->GetOutputPort(0));
  VTK_CREATE(vtkActor, actor);
  actor->SetMapper(mapper);
  VTK_CREATE(vtkRenderer, renderer);
  renderer->AddActor(actor);
  VTK_CREATE(vtkRenderWindow, renwin);
  renwin->AddRenderer(renderer);
  VTK_CREATE(vtkRenderWindowInteractor, iren);
  iren->SetRenderWindow(renwin);
  renwin->Render();
  iren->Start();
#endif

  CheckWinding(clipper);

  VTK_CREATE(vtkDataSetSurfaceFilter, surface);
  surface->SetInputConnection(clipper->GetOutputPort());
  surface->Update();

  if (surface->GetOutput()->GetNumberOfCells() != numSurfacePolys)
    {
    std::cout << "Expected " << numSurfacePolys << " triangles on the surface, got "
         << surface->GetOutput()->GetNumberOfCells() << std::endl;
    throw BoxClipTriangulateFailed();
    }
}

//-----------------------------------------------------------------------------

int BoxClipTriangulate(int, char *[])
{
  long seed = time(NULL);
  std::cout << "Random seed = " << seed << std::endl;
  vtkMath::RandomSeed(seed);
  vtkMath::Random();
  vtkMath::Random();
  vtkMath::Random();

  try
    {
    std::cout << "Checking triangle strip." << std::endl;
    Check2DPrimitive(VTK_TRIANGLE_STRIP, NumTriStripCells, TriStripCells);

    std::cout << "Checking quadrilaterals." << std::endl;
    Check2DPrimitive(VTK_QUAD, NumQuadCells, QuadCells);

    std::cout << "Checking pixels." << std::endl;
    Check2DPrimitive(VTK_PIXEL, NumPixelCells, PixelCells);

    std::cout << "Checking polygons." << std::endl;
    Check2DPrimitive(VTK_POLYGON, NumPolyCells, PolyCells);

    std::cout << "Checking hexahedrons." << std::endl;
    Check3DPrimitive(VTK_HEXAHEDRON, NumHexCells, HexCells,
                     NumExpectedHexSurfacePolys);

    std::cout << "Checking voxels." << std::endl;
    Check3DPrimitive(VTK_VOXEL, NumVoxelCells, VoxelCells,
                     NumExpectedVoxelSurfacePolys);

    std::cout << "Checking wedges." << std::endl;
    Check3DPrimitive(VTK_WEDGE, NumWedgeCells, WedgeCells,
                     NumExpectedWedgeSurfacePolys);

    std::cout << "Checking pyramids." << std::endl;
    Check3DPrimitive(VTK_PYRAMID, NumPyramidCells, PyramidCells,
                     NumExpectedPyramidSurfacePolys);
    }
  catch (BoxClipTriangulateFailed)
    {
    return EXIT_FAILURE;
    }

  return EXIT_SUCCESS;
}