diff --git a/CMake/vtkLegacyData.cmake b/CMake/vtkLegacyData.cmake
index 0921fced275bd3f56cc917918285a1da8bc0f0f8..865e91e3da80862b20efc5561f4e497251d89c66 100644
--- a/CMake/vtkLegacyData.cmake
+++ b/CMake/vtkLegacyData.cmake
@@ -3,6 +3,7 @@
# TODO: Reference testing data from each module only as needed.
set(data "DATA{${VTK_TEST_INPUT_DIR}/,REGEX:.*}")
foreach(d
+ FiberSurface
Infovis
Infovis/SQLite
Infovis/XML
diff --git a/Filters/Topology/CMakeLists.txt b/Filters/Topology/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8af0cc798165dfeb7edd245c616561b3e7095008
--- /dev/null
+++ b/Filters/Topology/CMakeLists.txt
@@ -0,0 +1,5 @@
+set(Module_SRCS
+ vtkFiberSurface.cxx
+ )
+
+vtk_module_library(vtkFiltersTopology ${Module_SRCS})
diff --git a/Filters/Topology/Testing/Cxx/CMakeLists.txt b/Filters/Topology/Testing/Cxx/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5e542b6de11047c8dc40137b20941fab6db076a9
--- /dev/null
+++ b/Filters/Topology/Testing/Cxx/CMakeLists.txt
@@ -0,0 +1,4 @@
+vtk_add_test_cxx(${vtk-module}CxxTests tests
+ TestFiberSurface.cxx,NO_VALID
+ )
+vtk_test_cxx_executable(${vtk-module}CxxTests tests)
diff --git a/Filters/Topology/Testing/Cxx/TestFiberSurface.cxx b/Filters/Topology/Testing/Cxx/TestFiberSurface.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..35ab1990e7f39a6a9c16527300b47a74c6a4867a
--- /dev/null
+++ b/Filters/Topology/Testing/Cxx/TestFiberSurface.cxx
@@ -0,0 +1,214 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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.
+
+=========================================================================*/
+
+
+// PURPOSE: FiberSurface test cases output generator (new vtkFiberSurface filter)
+
+#include <cstring>
+#include <sstream>
+#include <string>
+
+#include "vtkFiberSurface.h"
+#include "vtkTestUtilities.h"
+#include "vtkActor.h"
+#include "vtkCellArray.h"
+#include "vtkCleanPolyData.h"
+#include "vtkFloatArray.h"
+#include "vtkPolyData.h"
+#include "vtkPolyData.h"
+#include "vtkPolyDataReader.h"
+#include "vtkPolyDataWriter.h"
+#include "vtkSmartPointer.h"
+#include "vtkUnstructuredGrid.h"
+#include "vtkUnstructuredGridReader.h"
+
+int TestFiberSurface(int argc, char* argv[])
+{
+ bool pass = true;
+
+ /************** Input File paths using vtkTestUtilities ****************/
+ const char* inputDataFiles[3] = { vtkTestUtilities::ExpandDataFileName(
+ argc, argv, "Data/FiberSurface/one_cube.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv,
+ "Data/FiberSurface/one_cube_both_forking.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv,
+ "Data/FiberSurface/one_cube_closed.vtk") };
+
+ /************** FSCP File paths using vtkTestUtilities ****************/
+ const char* inputFSCPFiles[15] = { vtkTestUtilities::ExpandDataFileName(
+ argc, argv, "Data/FiberSurface/line_01.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_02.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_03.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_04.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_05.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_11.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_12.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_13.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_14.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_15.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_21.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_22.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_23.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_24.vtk"),
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/FiberSurface/line_25.vtk") };
+ /**********************************************/
+
+ /********************* Desired output arrays *********************/
+ std::string dataToCompare[15] = { "0.779,0.000,0.000,0.659,0.000,0.341,1.000,0.000,0.624,1.000,0."
+ "797,0.000,0.659,0.659,1.000,1.000,1.000,0.203,1.000,0.376,1."
+ "000,0.779,1.000,1.000",
+ "0.775,0.000,0.775,0.798,0.000,1.000,0.768,0.232,0.768,0.874,0.000,0.000,0.889,0.000,0.111,0."
+ "918,0.918,0.082,0.889,0.889,1.000,0.874,1.000,1.000,0.775,0.225,1.000",
+ "0.000,0.800,0.000,0.000,0.000,0.667,0.286,0.000,0.286,0.200,0.000,0.000,0.000,0.333,1.000,0."
+ "200,1.000,1.000,0.286,0.714,1.000,0.000,1.000,0.200",
+ "0.331,0.000,0.331,0.167,0.833,0.167,0.202,0.000,1.000,0.285,0.000,0.000,0.169,0.000,0.831,0."
+ "210,0.210,0.790,0.169,0.169,1.000,0.285,1.000,1.000,0.331,0.669,1.000",
+ "0.650,0.000,1.000,0.500,0.500,0.500,0.000,0.053,0.000,0.000,0.000,0.029,1.000,0.000,0.478,0."
+ "057,0.000,0.000,1.000,0.868,0.000,0.931,0.931,0.069,1.000,0.522,1.000,1.000,1.000,0.132,0.000,"
+ "0.971,1.000,0.057,1.000,1.000,0.000,1.000,0.947",
+ "0.500,0.500,1.000,1.000,1.000,0.714,1.000,0.667,1.000,0.714,1.000,1.000",
+ "0.000,0.901,0.000,0.233,0.767,0.233,0.000,0.271,0.729,0.886,0.000,0.000,0.241,0.000,0.759,1."
+ "000,0.000,0.263,1.000,0.201,0.000,0.064,0.064,1.000,1.000,1.000,0.460,1.000,0.376,1.000,0.469,"
+ "1.000,1.000,0.000,1.000,0.172",
+ "0.000,0.571,0.000,0.000,0.000,0.667,1.000,0.000,1.000,0.571,0.000,0.000,0.667,0.667,0.333,1."
+ "000,0.000,1.000,1.000,0.750,0.000,1.000,1.000,0.143,0.000,0.667,1.000,0.143,1.000,1.000,1.000,"
+ "0.000,1.000,0.000,1.000,0.750",
+ "0.000,0.250,0.000,0.000,0.000,0.167,0.250,0.000,0.250,0.100,0.000,0.000,0.000,0.833,1.000,0."
+ "100,1.000,1.000,0.250,0.750,1.000,0.000,1.000,0.750",
+ "0.333,0.000,0.000,0.333,0.000,0.667,1.000,0.000,0.667,0.333,0.333,0.667,1.000,0.333,0.667,1."
+ "000,1.000,0.667",
+ "0.000,0.000,0.300,0.300,0.000,0.300,0.300,0.700,0.300,0.000,0.700,0.300,0.700,0.000,0.300,1."
+ "000,0.000,0.300,0.700,0.700,0.300,1.000,0.700,0.300,1.000,0.700,1.000,0.700,0.700,1.000,0.300,"
+ "0.700,1.000,0.000,0.700,1.000",
+ "0.800,0.200,0.800,0.800,0.000,0.800,0.000,0.000,0.800,0.000,0.200,0.800,0.200,0.000,0.800,1."
+ "000,0.000,0.800,1.000,0.200,0.800,0.200,0.200,0.800,0.200,0.200,1.000,1.000,0.200,1.000,0.800,"
+ "0.200,1.000,0.000,0.200,1.000",
+ "0.828,0.000,0.828,0.737,0.000,1.000,0.856,0.144,0.856,1.000,0.000,0.828,1.000,0.144,0.856,1."
+ "000,0.172,1.000,0.856,0.144,1.000",
+ "0.000,0.739,0.000,0.000,0.000,0.146,0.427,0.000,0.427,0.854,0.000,0.146,1.000,0.000,0.427,1."
+ "000,0.739,0.000,0.854,0.854,1.000,1.000,1.000,0.261,1.000,0.573,1.000,0.261,1.000,1.000",
+ "0.977,0.023,0.977,0.980,0.000,0.980,0.000,0.000,0.671,0.000,0.363,0.637,0.329,0.000,0.671,1."
+ "000,0.000,0.980,1.000,0.023,0.977,0.363,0.363,0.637,0.329,0.329,1.000,1.000,0.020,1.000,0.977,"
+ "0.023,1.000,0.000,0.363,1.000" };
+
+ std::string outputString;
+ char firstFieldName[] = "f1";
+ char secondFieldName[] = "f2";
+ int cmpIndex = 0;
+
+ /************** output to error file ****************/
+ /****************************************************/
+
+ printf("FiberSurface test cases");
+ /****************************************************/
+
+ for (int i = 0; i < 3; i++) // loop for three input files
+ {
+
+ for (int j = 0; j < 5; j++) // loop for fifteen fscp files
+ {
+ // read and load .vtk input data file
+ vtkSmartPointer<vtkUnstructuredGridReader> mesh_reader =
+ vtkSmartPointer<vtkUnstructuredGridReader>::New();
+ mesh_reader->SetFileName(inputDataFiles[i]);
+ mesh_reader->Update();
+
+ // read and load .vtu file containing control polyline
+ vtkSmartPointer<vtkPolyDataReader> polyline_reader =
+ vtkSmartPointer<vtkPolyDataReader>::New();
+ polyline_reader->SetFileName(inputFSCPFiles[cmpIndex]);
+ polyline_reader->Update();
+
+ // extract polydata surface
+ vtkSmartPointer<vtkFiberSurface> fs = vtkSmartPointer<vtkFiberSurface>::New();
+ fs->SetInputData(0, mesh_reader->GetOutput());
+ fs->SetInputData(1, polyline_reader->GetOutput());
+ fs->SetField1(firstFieldName);
+ fs->SetField2(secondFieldName);
+ fs->Update();
+
+ /********************* cleaning and comparing FS Soup to desired output ********************/
+
+ vtkPolyData* polyData = fs->GetOutput();
+
+ /**** Duplicate Ponts ****/
+ int numOfPoints = polyData->GetNumberOfPoints();
+ int index = 0;
+ double* coords;
+
+ /**** NON Duplicate Ponts ****/
+ vtkSmartPointer<vtkCleanPolyData> clean = vtkSmartPointer<vtkCleanPolyData>::New();
+ clean->SetInputData(polyData);
+ clean->Update();
+ numOfPoints = clean->GetOutput()->GetNumberOfPoints();
+ std::cout << ".";
+
+ index = 0;
+ outputString = "";
+ std::stringstream ss;
+
+ // out non duplicate points
+ while (index < numOfPoints)
+ {
+ coords = clean->GetOutput()->GetPoint(index);
+ ss << std::fixed << std::setprecision(3) << coords[0] << "," << std::fixed
+ << std::setprecision(3) << coords[1] << "," << std::fixed << std::setprecision(3)
+ << coords[2];
+ if (index > 0)
+ outputString.append(",");
+ outputString.append(ss.str());
+ ss.str(std::string());
+ ss.clear();
+ index++;
+ } // end while
+
+ /**********************************************/
+ /*********** comparing the results ************/
+
+ if (strcmp(outputString.c_str(), dataToCompare[cmpIndex].c_str()) != 0)
+ {
+ pass = false;
+ { // writing to file
+ printf("\n\n/**************************************/");
+ printf("\n/********Test Unsuccessful*************/");
+ printf("\nInput Data: %s", (std::string(inputDataFiles[i]).substr(6, 31)).c_str());
+ printf(
+ "\nInput FSCP: %s", (std::string(inputFSCPFiles[cmpIndex]).substr(6, 30)).c_str());
+ printf("\nString to compare: %s", dataToCompare[cmpIndex].c_str());
+ printf("\nOutput String : %s", outputString.c_str());
+ printf("\n/**************************************/");
+ } // end of writing
+ }
+ cmpIndex++;
+ } // end for
+ } // end for
+
+ for (int i = 0; i < 15; ++i)
+ {
+ delete [] inputFSCPFiles[i];
+ }
+ delete [] inputDataFiles[0];
+ delete [] inputDataFiles[1];
+ delete [] inputDataFiles[2];
+
+ if (pass)
+ {
+ cout << "\nTest Successful!!!" << endl;
+ return EXIT_SUCCESS;
+ }
+
+ cout << "\nTest Unsuccessful." << endl;
+ return EXIT_FAILURE;
+}
diff --git a/Filters/Topology/module.cmake b/Filters/Topology/module.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..a47982435c277b2e46dd998bb2132ebcc0dd836c
--- /dev/null
+++ b/Filters/Topology/module.cmake
@@ -0,0 +1,13 @@
+vtk_module(vtkFiltersTopology
+ GROUPS
+ StandAlone
+ TEST_DEPENDS
+ vtkTestingRendering
+ vtkIOLegacy
+ KIT
+ vtkFilters
+ DEPENDS
+ vtkCommonCore
+ vtkCommonDataModel
+ vtkCommonExecutionModel
+ )
diff --git a/Filters/Topology/vtkFiberSurface.cxx b/Filters/Topology/vtkFiberSurface.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2265e1e667947ecad2683dfd82d32ffcd8103ae8
--- /dev/null
+++ b/Filters/Topology/vtkFiberSurface.cxx
@@ -0,0 +1,1075 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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 "vtkFiberSurface.h"
+
+#include "vtkCell.h"
+#include "vtkCellArray.h"
+#include "vtkDoubleArray.h"
+#include "vtkInformation.h"
+#include "vtkInformationVector.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkUnstructuredGrid.h"
+
+// lookup table for powers of 3 shifts in the marching tetrahedra cases
+// which is described in greyTetTriangles table.
+// remind that we use 0, 1 and 2 to represent (W)hite, (G)rey and (B)lack cases. For each
+// tetrahedron, cases for four vertices can be represented
+// by a four-digit number, such as 0001. We assume that all vertices are in CCW order.
+// The array greyTetTriangles actually records all of 81 such cases. The order of case
+// index starts from right-most to the left-most digit: starting from 0000 to 0002,
+// then from 0010 to 0022, then from 0100 to 0222, finally from 1000 to 2222.
+// It is easy to observe that:
+// from 0001 to 0002, as case number in the first dight is incremented by 1,
+// we only need to skip 1 index in the greyTetTriangles table.
+//
+// from 0010 to 0020, as case number in the second dight is incremented by 1,
+// we need to skip 3 indices (0010, 0011, 0012, 0020)
+//
+// from 0100 to 0200, as case number in the third dight is incremented by 1,
+// we need to skip 9 indices (0100, 0101, 0102, 0110, 0111,
+// 0112, 0120, 0121, 0122, 0200)
+//
+// from 1000 to 2000, as case number in the fourth dight is incremented by 1,
+// we need to skip 27 indices (1000, 1001, 1002, 1010, 1011,
+// 1012, 1020, 1021, 1022, 1100,
+// 1101, 1102, 1110, 1111, 1112,
+// 1120, 1121, 1122, 1200, 1201,
+// 1202, 1210, 1211, 1212, 1220,
+// 1221, 1222, 2000)
+// Given case classifications for four vertices in a tetrahedron, this ternaryShift array
+// could be used to quickly locate the index number in the marching tetrahedron case
+// table. This array can also be used in the clipping case look-up table
+// clipTriangleVertices.
+static const char ternaryShift[4] = { 1, 3, 9, 27 };
+
+//----------------------------------------------------------------------------
+
+// In the Marching Tethrhedron with Grey case, the iso-surface can be either a triangle
+// ,quad or null. The number of triangles in each case is at most 2. This array
+// records the number of triangles for every case.
+static const int nTriangles[81] = {
+ 0, 0, 1, 0, 0, 1, 1, 1, 2, // cases 0000-0022
+ 0, 0, 1, 0, 1, 1, 1, 1, 1, // cases 0100-0122
+ 1, 1, 2, 1, 1, 1, 2, 1, 1, // cases 0200-0222
+
+ 0, 0, 1, 0, 1, 1, 1, 1, 1, // cases 1000-1022
+ 0, 1, 1, 1, 0, 1, 1, 1, 0, // cases 1100-1122
+ 1, 1, 1, 1, 1, 0, 1, 0, 0, // cases 1200-1222
+
+ 1, 1, 2, 1, 1, 1, 2, 1, 1, // cases 2000-2022
+ 1, 1, 1, 1, 1, 0, 1, 0, 0, // cases 2100-2122
+ 2, 1, 1, 1, 0, 0, 1, 0, 0 // cases 2200-2222
+};
+
+//----------------------------------------------------------------------------
+
+// array of vertices for triangles in the marching tetrahedron cases
+// each vertex on the tetra is marked as (B)lack, (W)hite or (G)rey
+// there are total 81 cases. Each case contains at most two triangles.
+// The order these cases arranged are as follows: starting from 0000 to 0002,
+// then from 0010 to 0022, then from 0100 to 0222, finally from 1000 to 2222.
+static const vtkFiberSurface::BaseVertexType greyTetTriangles[81][2][3] = {
+ { // 0. case 0000 (A)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 1. case 0001 (B)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 2. case 0002 (D)
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 3. case 0010 (B)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 4. case 0011 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 5. case 0012 (F)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 6. case 0020 (D)
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 7. case 0021 (F)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 8. case 0022 (I)
+ { vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_12 } },
+ { // 9. case 0100 (B)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 10. case 0101 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 11. case 0102 (F)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 12. case 0110 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 13. case 0111 (E)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_2 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 14. case 0112 (H)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 15. case 0120 (F)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 16. case 0121 (H)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 17. case 0122 (K)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 18. case 0200 (D)
+ { vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 19. case 0201 (F)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 20. case 0202 (I)
+ { vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_01 } },
+ { // 21. case 0210 (F)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 22. case 0211 (H)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 23. case 0212 (K)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 24. case 0220 (I)
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_02 } },
+ { // 25. case 0221 (K)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 26. case 0222 (M)
+ { vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 27. case 1000 (B)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 28. case 1001 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 29. case 1002 (F)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 30. case 1010 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 31. case 1011 (E)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_1 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 32. case 1012 (H)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 33. case 1020 (F)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 34. case 1021 (H)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 35. case 1022 (K)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 36. case 1100 (C)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 37. case 1101 (E)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_3 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 38. case 1102 (H)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 39. case 1110 (E)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_2 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 40. case 1111 (G)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 41. case 1112 (J)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_3 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 42. case 1120 (H)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 43. case 1121 (J)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_2 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 44. case 1122 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 45. case 1200 (F)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 46. case 1201 (H)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 47. case 1202 (K)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 48. case 1210 (H)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 49. case 1211 (J)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_3 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 50. case 1212 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 51. case 1220 (K)
+ { vtkFiberSurface::bv_vertex_3, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 52. case 1221 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 53. case 1222 (N)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 54. case 2000 (D)
+ { vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 55. case 2001 (F)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 56. case 2002 (I)
+ { vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_02 } },
+ { // 57. case 2010 (F)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 58. case 2011 (H)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 59. case 2012 (K)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_23 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 60. case 2020 (I)
+ { vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_23 } },
+ { // 61. case 2021 (K)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 62. case 2022 (M)
+ { vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_23, vtkFiberSurface::bv_edge_12 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 63. case 2100 (F)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 64. case 2101 (H)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 65. case 2102 (K)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_01 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 66. case 2110 (H)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 67. case 2111 (J)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_vertex_1 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 68. case 2112 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 69. case 2120 (K)
+ { vtkFiberSurface::bv_vertex_2, vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_03 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 70. case 2121 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 71. case 2122 (N)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 72. case 2200 (I)
+ { vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_edge_13, vtkFiberSurface::bv_edge_02, vtkFiberSurface::bv_edge_12 } },
+ { // 73. case 2201 (K)
+ { vtkFiberSurface::bv_vertex_0, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 74. case 2202 (M)
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_12, vtkFiberSurface::bv_edge_13 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 75. case 2210 (K)
+ { vtkFiberSurface::bv_vertex_1, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 76. case 2211 (L)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 77. case 2212 (N)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 78. case 2220 (M)
+ { vtkFiberSurface::bv_edge_01, vtkFiberSurface::bv_edge_03, vtkFiberSurface::bv_edge_02 },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 79. case 2221 (N)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } },
+ { // 80. case 2222 (O)
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used },
+
+ { vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used, vtkFiberSurface::bv_not_used } }
+};
+
+//----------------------------------------------------------------------------
+
+// conversion from the enum semantics for edges to actual edge numbers
+// depends on the ordering of bv_edge_** in BaseVertexType enum
+static const int edge2endpoints[6][2] = { { 0, 1 }, { 0, 2 }, { 0, 3 }, { 1, 2 }, { 1, 3 },
+ { 2, 3 } };
+
+//----------------------------------------------------------------------------
+
+// convert edge_*_param_* enum to edge numbers
+// depends on the ordering of edge_0 and edge_1 enums (i.e. edge_0 + 2 == edge_1 + 1
+// == edge_2)
+static const int clip2points[3][2] = { { 1, 2 }, { 2, 0 }, { 0, 1 } };
+
+//----------------------------------------------------------------------------
+
+// this table lists the number of triangles per case for fiber clipping
+static const int nClipTriangles[27] = {
+ 0, 1, 2, 1, 2, 3, 2, 3, 2, // cases 000 - 022
+ 1, 2, 3, 2, 1, 2, 3, 2, 1, // cases 100 - 122
+ 2, 3, 2, 3, 2, 1, 2, 1, 0 // cases 200 - 222
+};
+
+//----------------------------------------------------------------------------
+
+// with up to three triangles, we can have up to 9 vertices specified
+// note that this may lead to redundant interpolation (as in MC/MT), but we gain in
+// clarity by doing it this way.
+// This array therefore specifies the vertices of each triangle to be rendered in the
+// clipping process.
+static const int clipTriangleVertices[27][3][3] = {
+ // 0. case 000: A - empty
+ { { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 1. case 001: B - point-triangle
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 2. case 002: D - stripe
+ { { vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_1_parm_0,
+ vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_1_parm_1,
+ vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 3. case 010: B - point-triangle
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 4. case 011: C - edge-quad
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 5. case 012: E - point-stripe
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_0_parm_0,
+ vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_0_parm_0,
+ vtkFiberSurface::edge_1_parm_0 } },
+ // 6. case 020: D - stripe
+ { { vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_2_parm_0,
+ vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_2_parm_1,
+ vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 7. case 021: E - point-stripe
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_2_parm_1,
+ vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_2_parm_1,
+ vtkFiberSurface::edge_0_parm_1 } },
+ // 8. case 022: D - stripe
+ { { vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_0_parm_1,
+ vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_0_parm_0,
+ vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 9. case 100: B - point-triangle
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 10. case 101: C - edge-quad
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 11. case 102: E - point-stripe
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::edge_0_parm_0, vtkFiberSurface::edge_1_parm_1,
+ vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_1_parm_1,
+ vtkFiberSurface::edge_2_parm_1 } },
+ // 12. case 110: C - edge-quad
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 13. case 111: F - entire triangle
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::vertex_1, vtkFiberSurface::vertex_2 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 14. case 112: C - edge-quad
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 15. case 120: E - point-stripe
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_1_parm_0,
+ vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_1_parm_0,
+ vtkFiberSurface::edge_2_parm_0 } },
+ // 16. case 121: C - edge-quad
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 17. case 122: B - point-triangle
+ { { vtkFiberSurface::vertex_2, vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 18. case 200: D - stripe
+ { { vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_0_parm_0,
+ vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_0_parm_1,
+ vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 19. case 201: E - point-stripe
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::edge_1_parm_1, vtkFiberSurface::edge_2_parm_0,
+ vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_2_parm_0,
+ vtkFiberSurface::edge_0_parm_0 } },
+ // 20. case 202: D - stripe
+ { { vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_2_parm_1,
+ vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_2_parm_0,
+ vtkFiberSurface::edge_0_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 21. case 210: E - point-stripe
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::edge_2_parm_0, vtkFiberSurface::edge_0_parm_1,
+ vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::edge_1_parm_0, vtkFiberSurface::edge_0_parm_1,
+ vtkFiberSurface::edge_1_parm_1 } },
+ // 22. case 211: C - edge-quad
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_1 },
+
+ { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 23. case 212: B - point-triangle
+ { { vtkFiberSurface::vertex_1, vtkFiberSurface::edge_0_parm_1, vtkFiberSurface::edge_2_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 24. case 220: D - stripe
+ { { vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_1_parm_1,
+ vtkFiberSurface::edge_1_parm_0 },
+
+ { vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_1_parm_0,
+ vtkFiberSurface::edge_2_parm_0 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 25. case 221: B - point-triangle
+ { { vtkFiberSurface::vertex_0, vtkFiberSurface::edge_2_parm_1, vtkFiberSurface::edge_1_parm_1 },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } },
+ // 26. case 222: A - empty
+ { { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used },
+
+ { vtkFiberSurface::not_used, vtkFiberSurface::not_used, vtkFiberSurface::not_used } }
+};
+
+//----------------------------------------------------------------------------
+
+vtkStandardNewMacro(vtkFiberSurface);
+
+//----------------------------------------------------------------------------
+
+vtkFiberSurface::vtkFiberSurface()
+{
+ // number of input ports is 2
+ this->SetNumberOfInputPorts(2);
+ this->Fields[0] = this->Fields[1] = NULL;
+}
+
+//----------------------------------------------------------------------------
+
+void vtkFiberSurface::SetField1(char* nm)
+{
+ this->Fields[0] = nm;
+}
+
+//----------------------------------------------------------------------------
+
+void vtkFiberSurface::SetField2(char* nm)
+{
+ this->Fields[1] = nm;
+}
+
+//----------------------------------------------------------------------------
+
+void vtkFiberSurface::PrintSelf(ostream& os, vtkIndent indent)
+{
+ this->Superclass::PrintSelf(os, indent);
+}
+
+//----------------------------------------------------------------------------
+
+int vtkFiberSurface::FillInputPortInformation(int port, vtkInformation* info)
+{
+ switch (port)
+ {
+ // port 0 expects a tetrahedral mesh as input data
+ case 0:
+ info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
+ return 1;
+ // port 1 expects a fiber surface control polygon (FSCP)
+ case 1:
+ info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkPolyData");
+ return 1;
+ }
+ return 0;
+}
+
+//----------------------------------------------------------------------------
+
+int vtkFiberSurface::RequestData(vtkInformation* vtkNotUsed(request),
+ vtkInformationVector** inputVector, vtkInformationVector* outputVector)
+{
+ // obtain the input/output port info
+ vtkInformation* inMeshInfo = inputVector[0]->GetInformationObject(0);
+ vtkInformation* inLinesInfo = inputVector[1]->GetInformationObject(0);
+ vtkInformation* outInfo = outputVector->GetInformationObject(0);
+
+ // get the input regular grid and fiber surface control polygon (FSCP)
+ vtkUnstructuredGrid* mesh =
+ vtkUnstructuredGrid::SafeDownCast(inMeshInfo->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData* lines = vtkPolyData::SafeDownCast(inLinesInfo->Get(vtkDataObject::DATA_OBJECT()));
+ vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
+
+ // get the dataset statistics
+ vtkIdType numCells = mesh->GetNumberOfCells();
+ vtkIdType numPts = mesh->GetNumberOfPoints();
+ vtkIdType numArrays = mesh->GetPointData()->GetNumberOfArrays();
+
+ // check if the data set is empty or not.
+ // check if the data set contains at least two scalar arrays.
+ if (numCells < 1 || numPts < 1 || numArrays < 2)
+ {
+ vtkErrorMacro("No input data. Two fields are required for fiber surface generation");
+ return 1;
+ }
+
+ // check if two scalar fields are specified by the user.
+ // if not specified, an error message will be shown.
+ if (!this->Fields[0] || !this->Fields[1])
+ {
+ vtkErrorMacro("Two scalar fields need to be specified.");
+ return 1;
+ }
+ else
+ {
+ vtkSmartPointer<vtkDataArray> array;
+ for (int i = 0; i < 2; i++)
+ {
+ array = mesh->GetPointData()->GetArray(this->Fields[i]);
+ if (!array)
+ {
+ vtkErrorMacro("Names of the scalar array do not exist.");
+ return 1;
+ }
+ }
+ }
+
+ // extract the points of the subdivided tetrahedra
+ vtkPoints* meshPoints = mesh->GetPoints();
+
+ // allocate points and cell storage for computing fiber surface structure
+ vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
+ vtkSmartPointer<vtkCellArray> newPolys = vtkSmartPointer<vtkCellArray>::New();
+
+ // extract two scalar field arrays and put into one structure
+ vtkDataArray* fieldScalars[2] = { mesh->GetPointData()->GetArray(this->Fields[0]),
+ mesh->GetPointData()->GetArray(this->Fields[1]) };
+
+ // extract a fiber surface for every edge in the FSCP.
+ // if the FSCP has no edges, this loop will not start.
+ // Algorithm:
+ // 1. Extract base fiber surface using marching tetrahedron
+ // 2. Clip the base fiber surface to extract the exact fiber surface with respect to
+ // each line segment in FSCP
+ const vtkIdType numberOfLines = lines->GetNumberOfCells();
+ for (vtkIdType lineIndex = 0; lineIndex < numberOfLines; ++lineIndex)
+ {
+ // For each line segment of the FSCP
+ vtkCell* line = lines->GetCell(lineIndex);
+
+ // test if the current cell of FSCP is a line or not.
+ // The computation only proceed if the current cell is a line.
+ if (line->GetNumberOfPoints() != 2)
+ {
+ vtkWarningMacro("Current cell index " << lineIndex << " in the FSCP is not of a line type.");
+ continue;
+ }
+
+ // get the starting point of the line
+ double pointStart[3], pointEnd[3];
+ lines->GetPoints()->GetPoint(line->GetPointId(0), pointStart);
+
+ // get the end point of the line
+ lines->GetPoints()->GetPoint(line->GetPointId(1), pointEnd);
+
+ // first point is the origin of the parametric form of line
+ const double origin[2] = { pointStart[0], pointStart[1] };
+
+ // if the input point coordinates are normalised values. We need to interpolate these
+ // values to the actual scalar values.
+ // obtain the direction vector of the line segment
+ const double direction[2] = { pointEnd[0] - origin[0], pointEnd[1] - origin[1] };
+
+ // compute length of the line segment
+ const double length = sqrt(direction[0] * direction[0] + direction[1] * direction[1]);
+
+ // if the length of the current line is zero, then skip to the next cell.
+ if (length == 0.0f)
+ {
+ vtkWarningMacro("End points of the current line index "
+ << lineIndex << " in the FSCP colocate on the same point.");
+ continue;
+ }
+
+ // compute the normal vector to the line segment
+ const double normal[2] = { direction[1] / length, -direction[0] / length };
+
+ // given a line segment with one of its endpoint origin and its normal vector normal
+ // given an arbitrary point p
+ // the signed distance from p to the line can be computed using the Hess Normal Form:
+ // signedDistance = dot(p - origin, normal)
+ // = dot(p,normal) - dot(origin,normal)
+ // since dot(origin,normal) is an invariant, therefore we compute it first to avoid
+ // duplicate computation in the following steps.
+ const double dotOriginNormal = normal[0] * origin[0] + normal[1] * origin[1];
+
+ // get the number of tetras in the domain
+ const vtkIdType numberOfTets = mesh->GetNumberOfCells();
+
+ // iterate through every cell of the domain and extract its fiber surface.
+ // note that each cell is a tetrahedron.
+ for (vtkIdType tetIndex = 0; tetIndex < numberOfTets; ++tetIndex)
+ {
+ // update progress of the extraction
+ this->UpdateProgress((tetIndex + 1.0) / numberOfTets);
+
+ // obtain the current tetra cell
+ vtkCell* tet = mesh->GetCell(tetIndex);
+
+ // check if the current cell is a tetrahedron type or not
+ // if not, skip this cell.
+ if (mesh->GetCellType(tetIndex) != VTK_TETRA || tet->GetNumberOfPoints() != 4)
+ {
+ vtkWarningMacro("Current cell " << tetIndex << "is not of a tetrahedron type.");
+ continue;
+ }
+
+ // case number for the current tetra cell, initially set to zero
+ unsigned caseNumber = 0;
+
+ // classify four vertices of the tetra with respect to the signed distance to the
+ // line
+ double distancesToLine[4];
+
+ // for each vertex in the tetra cell
+ for (int vertexIndex = 0; vertexIndex < 4; ++vertexIndex)
+ {
+ // get the Id of the vertex of the tetra
+ const vtkIdType pointId = tet->GetPointId(vertexIndex);
+
+ // compute signed distance between the image of tetra vertex in the range
+ // and the control line using Hessian Normal Form:
+ // signedDistance = dot(p - origin, normal)
+ // = dot(p,normal) - dot(origin,normal)
+ distancesToLine[vertexIndex] = fieldScalars[0]->GetTuple1(pointId) * normal[0] +
+ fieldScalars[1]->GetTuple1(pointId) * normal[1] - dotOriginNormal;
+
+ // classify the tetra vertex based on the sign of the distance
+ // if distance == 0 : then p is on the line
+ // if distance > 0 : then p is on the right side
+ // if distance < 0; then p is on the left side
+ if (distancesToLine[vertexIndex] == 0.0)
+ {
+ // locate the index number using ternaryShift array.
+ // please see the comment on this array.
+ caseNumber += ternaryShift[vertexIndex];
+ }
+ else if (distancesToLine[vertexIndex] > 0.0)
+ {
+ // locate the index number using ternaryShift array.
+ // please see the comment on this array.
+ caseNumber += 2 * ternaryShift[vertexIndex];
+ }
+ }
+
+ // extract the base fiber surface using Marching Tetrahedra
+ // loop starts only when there is at least one triangle in this case.
+ for (int triangleIndex = 0; triangleIndex < nTriangles[caseNumber]; ++triangleIndex)
+ {
+ // Coordinates for each triangle points
+ double trianglePoints[3][3];
+
+ // Clipping parameter for the base fiber surface
+ double triangleParameters[3];
+
+ // Clipping case number is initially set to zero.
+ unsigned triangleCaseNumber = 0;
+
+ // global index in the input data of the current vertex
+ vtkIdType dataIndex;
+
+ // scalar values of the two fields of the current vertex;
+ double pointField1, pointField2;
+
+ // For each vertex in the base fiber surface
+ for (int pointIndex = 0; pointIndex < 3; ++pointIndex)
+ {
+ // get the type of the triangle vertex from the lookup table
+ const int type = greyTetTriangles[caseNumber][triangleIndex][pointIndex];
+ switch (type)
+ {
+ case bv_vertex_0:
+ case bv_vertex_1:
+ case bv_vertex_2:
+ case bv_vertex_3:
+ // if the triangle vertex coincides with the tetra vertex,
+ // there will be a grey case.
+ // copy grey vertex to output triangle point array
+ double point[3];
+ dataIndex = tet->GetPointId(type - bv_vertex_0);
+ meshPoints->GetPoint(dataIndex, point);
+ memcpy(trianglePoints[pointIndex], point, sizeof(point));
+
+ // get the current vertex scalar values
+ pointField1 = fieldScalars[0]->GetTuple1(dataIndex);
+ pointField2 = fieldScalars[1]->GetTuple1(dataIndex);
+
+ // compute parameter on the parametric line for each triangle point.
+ // Assume edgeRange is a vector to represent the vector beteen interpolated
+ // range values and origin of polygon line edge.
+ // Assume direction is the direction vector of the current polygon edge.
+ // The projection of the range values onto the polygon edge
+ // can be computed as:
+ // Proj^direction_edgeRange = dot(edgeRange,direction) / |direction|^2
+ // this projection is the parameter t used in the clipping case.
+ // if t < 0 or t > 1: then the current vertex is outside the current line
+ // segment of FSCP
+ // if 0 << t << 1 : then the current vertex is within the current line
+ // segment of FSCP
+ triangleParameters[pointIndex] = ((pointField1 - origin[0]) * direction[0] +
+ (pointField2 - origin[1]) * direction[1]) /
+ (length * length);
+
+ if (triangleParameters[pointIndex] > 1.0)
+ {
+ // locate the index number in the clipping case table
+ triangleCaseNumber += 2 * ternaryShift[pointIndex];
+ }
+ else if (triangleParameters[pointIndex] >= 0.0)
+ {
+ // locate the index number in the clipping case table
+ triangleCaseNumber += ternaryShift[pointIndex];
+ }
+ // less than zero case adds nothing to triangleCaseNumber
+ break;
+
+ case bv_edge_01:
+ case bv_edge_02:
+ case bv_edge_03:
+ case bv_edge_12:
+ case bv_edge_13:
+ case bv_edge_23:
+ {
+ // If the triangle vertex happens on the tetra edges,
+ // compute interpolation mixing value.
+ // Note that (type - bv_edge_01) represents the current edge index.
+ // Assume an edge with two end points u and v.
+ // Since we are now in edge case, given the signed distances of u and v,
+ // the alpha value can be computed as:
+ // alpha = signedDistance(u) / (signedDistance(u) - signedDistance(v))
+ const double alpha = distancesToLine[edge2endpoints[type - bv_edge_01][0]] /
+ (distancesToLine[edge2endpoints[type - bv_edge_01][0]] -
+ distancesToLine[edge2endpoints[type - bv_edge_01][1]]);
+
+ // convert enum to pair of endpoints and get their id in the point set
+ const vtkIdType pointIds[2] = { tet->GetPointId(edge2endpoints[type - bv_edge_01][0]),
+ tet->GetPointId(edge2endpoints[type - bv_edge_01][1]) };
+
+ // get coordinates of the edge end points
+ double point0[3], point1[3];
+ meshPoints->GetPoint(pointIds[0], point0);
+ meshPoints->GetPoint(pointIds[1], point1);
+
+ // compute triangle vertex coordinates using linear interpolation
+ trianglePoints[pointIndex][0] = (1.0 - alpha) * point0[0] + alpha * point1[0];
+ trianglePoints[pointIndex][1] = (1.0 - alpha) * point0[1] + alpha * point1[1];
+ trianglePoints[pointIndex][2] = (1.0 - alpha) * point0[2] + alpha * point1[2];
+
+ // compute triangle vertex range values using linear interpolation.
+ const double edgeFields[2] = { (1.0 - alpha) *
+ fieldScalars[0]->GetTuple1(pointIds[0]) +
+ alpha * fieldScalars[0]->GetTuple1(pointIds[1]),
+ (1.0 - alpha) * fieldScalars[1]->GetTuple1(pointIds[0]) +
+ alpha * fieldScalars[1]->GetTuple1(pointIds[1]) };
+
+ // compute parameter on the parametric line for each triangle point.
+ // Assume edgeRange is a vector to represent the vector beteen interpolated
+ // range values and origin of polygon line edge.
+ // Assume direction is the direction vector of the current polygon edge.
+ // The projection of the range values onto the polygon edge
+ // can be computed as:
+ // Proj^direction_edgeRange = dot(edgeRange,direction) / |direction|^2
+ // this projection is the parameter t used in the clipping case.
+ // if t < 0 or t > 1: then the current vertex is outside the current line
+ // segment of FSCP
+ // if 0 << t << 1 : then the current vertex is within the current line
+ // segment of FSCP
+ triangleParameters[pointIndex] = ((edgeFields[0] - origin[0]) * direction[0] +
+ (edgeFields[1] - origin[1]) * direction[1]) /
+ (length * length);
+
+ if (triangleParameters[pointIndex] > 1.0)
+ {
+ // locate the index number in the clipping case table
+ triangleCaseNumber += 2 * ternaryShift[pointIndex];
+ }
+ else if (triangleParameters[pointIndex] >= 0.0)
+ {
+ // locate the index number in the clipping case table
+ triangleCaseNumber += ternaryShift[pointIndex];
+ }
+ // less than zero case adds nothing to triangleCaseNumber
+ break;
+ }
+ default:
+ // report error in case an invalid triangle is being extracted
+ vtkErrorMacro(<< "Invalid value in the marching tetrahedra case: " << caseNumber);
+ break;
+ }
+ }
+ // clip or cull the triangle from the base fiber surface.
+ int counter = 0;
+ vtkIdType pts[3];
+ for (int tindex = 0; tindex != nClipTriangles[triangleCaseNumber];
+ ++tindex)
+ {
+ for (int vertexIndex = 0; vertexIndex != 3; ++vertexIndex)
+ {
+ // array to holds indices of the two edge endpoints
+ int ids[2];
+
+ // interpolant value of the triangle vertex with respect to the control line
+ double alpha;
+ const int type = clipTriangleVertices[triangleCaseNumber][tindex][vertexIndex];
+ switch (type)
+ {
+ case vertex_0:
+ case vertex_1:
+ case vertex_2:
+ pts[counter] = newPoints->InsertNextPoint(trianglePoints[type - vertex_0]);
+ break;
+ case edge_0_parm_0:
+ case edge_1_parm_0:
+ case edge_2_parm_0:
+ ids[0] = clip2points[type - edge_0_parm_0][0];
+ ids[1] = clip2points[type - edge_0_parm_0][1];
+
+ // compute which endpoint we are clipping and its interpolant
+ alpha = (0.0 - triangleParameters[ids[0]]) /
+ (triangleParameters[ids[1]] - triangleParameters[ids[0]]);
+
+ // interpolate point position to clipping corner
+ pts[counter] = newPoints->InsertNextPoint(
+ (1.0 - alpha) * trianglePoints[ids[0]][0] + alpha * trianglePoints[ids[1]][0],
+ (1.0 - alpha) * trianglePoints[ids[0]][1] + alpha * trianglePoints[ids[1]][1],
+ (1.0 - alpha) * trianglePoints[ids[0]][2] + alpha * trianglePoints[ids[1]][2]);
+ break;
+ case edge_0_parm_1:
+ case edge_1_parm_1:
+ case edge_2_parm_1:
+ ids[0] = clip2points[type - edge_0_parm_1][0];
+ ids[1] = clip2points[type - edge_0_parm_1][1];
+
+ // compute which endpoint we are clipping and its interpolant
+ alpha = (1.0 - triangleParameters[ids[0]]) /
+ (triangleParameters[ids[1]] - triangleParameters[ids[0]]);
+
+ // interpolate point position to clipping corner
+ pts[counter] = newPoints->InsertNextPoint(
+ (1.0 - alpha) * trianglePoints[ids[0]][0] + alpha * trianglePoints[ids[1]][0],
+ (1.0 - alpha) * trianglePoints[ids[0]][1] + alpha * trianglePoints[ids[1]][1],
+ (1.0 - alpha) * trianglePoints[ids[0]][2] + alpha * trianglePoints[ids[1]][2]);
+ break;
+ default:
+ vtkErrorMacro(<< "Invalid value in clipping triangle case: " << triangleCaseNumber);
+ }
+ if (++counter == 3)
+ {
+ newPolys->InsertNextCell(3, pts);
+ counter = 0;
+ }
+ }
+ }
+ }
+ }
+ }
+ // store the fiber surface structure to the output polyData
+ output->SetPoints(newPoints);
+ output->SetPolys(newPolys);
+ return 1;
+}
diff --git a/Filters/Topology/vtkFiberSurface.h b/Filters/Topology/vtkFiberSurface.h
new file mode 100644
index 0000000000000000000000000000000000000000..5635680e257e574733c93cf2ee04f335e6d779a8
--- /dev/null
+++ b/Filters/Topology/vtkFiberSurface.h
@@ -0,0 +1,400 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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.
+
+=========================================================================*/
+/**
+* @class vtkFiberSurface
+* @brief Given a fiber surface control polygon (FSCP) and an
+* unstructured grid composed of tetrahedral cells with two scalar arrays, this filter
+* computes the corresponding fiber surfaces.
+*
+* @section Introduction
+* Fiber surfaces are constructed from sets of fibers, the multivariate analogues
+* of isolines. The original paper [0] offers a general purpose method that produces
+* separating surfaces representing boundaries in bivariate fields. This filter is based
+* on an improvement over [0] which computes accurate and exact fiber surfaces. It can
+* handle arbitrary input polygons including open polygons or self-intersecting polygons.
+* The current implementation can better captures sharp features induced by polygon
+* bends [1].
+*
+* [0] Hamish Carr, Zhao Geng, Julien Tierny, Amit Chattopadhyay and Aaron Knoll,
+* Fiber Surfaces: Generalizing Isosurfaces to Bivariate Data,
+* Computer Graphics Forum, Volume 34, Issue 3, Pages 241-250, (EuroVis 2015)
+*
+* [1] Pavol Klacansky, Julien Tierny, Hamish Carr, Zhao Geng,
+* Fast and Exact Fiber Surfaces for Tetrahedral Meshes,
+* Paper in submission, 2015
+*
+* @section Algorithm For Extracting An Exact Fiber Surface
+* Require: R.1 A 3D domain space represented by an unstructured grid composed of
+* tetrahedral cells
+* R.2 Two scalar fields, f1 and f2, that map the domain space to a 2D range
+* space. These fields are assumed to be known at vertices of the
+* unstructured grid: i.e they are two fields associated with the
+* unstructured grid.
+* R.3 A Fiber Surface Control Polygon (FSCP) defined in the range space as a
+* list of line segments (l1, l2, ..., ln). The FSCP may be an open polyline
+* or a self-intersecting polygon.
+*
+* 1. For each line segment l in FSCP, we ignore the endpoints of the line and assume
+* this line extends to infinity. This line will then separate the range and its
+* inverse image, i.e fiber surfaces, will also separate the domain. Based on the
+* signed distance d between the image of a cell vertex v and line l in the range,
+* v can be classified as white (d < 0), grey (d == 0) or black (d>0). The
+* interpolation parameter between two vertices v1 and v2 in a cell edge can be
+* computed as |d1| / (|d2|+|d1|), where d1 and d2 are the signed distances between
+* images of v1,v2 and line l in the range. Once the classification and interpolation
+* parameters for all vertices in a cell are known, then we can apply the Marching
+* Tetrahedra algorithm. For each tetrahedron, this produces a planar cut which we
+* refer to as a base fiber surface.
+*
+* 2. After generating the base fiber surface in each cell, we need a further clipping
+* process to obtain the accurate fiber surface. Clipping is based on classifying the
+* vertices of each triangle as follows:
+* Given a line segment in the fiber surface control polygon (FSCP) parameterised
+* from 0 to 1, we know that every triangle vertex in the base fiber surface belongs
+* to the fiber surface, and that the fiber through each vertex can be parameterised
+* with respect to the line segment. Therefore, we compute the parameter t for each
+* vertex and use it to classify the vertex as:
+* 0: t < 0 Vertex is below the clipping range [0,1] and will be removed
+* 1: 0 ≤ t ≤ 1 Vertex is inside the clipping range [0,1] and is retained
+* 2: 1 < t Vertex is above the clipping range [0,1] and will be removed
+* Based on the classification, we can further clip the triangle to obtain the final
+* surface.
+*
+* 3. Repeating steps 1 and 2 for every line segment in FSCP and iterating through each
+* cell will generate the final fiber surfaces in the domain.
+*
+* @section VTK Filter Design
+* As stated in part B (R.1), we will compute the fiber surface over an unstructured grid.
+* This grid will have to be an input of the VTK filter. The two scalar fields, however,
+* are assumed to be stored in the VTK unstructured grid, and will be specified using the
+* SetField1() and SetField2() accessors. The FSCP will be passed into the filter as a
+* second input port. The data type of FSCP is required to be a vtkPolyData, with each of
+* its cell defined as a line segment and its point coordinates defined in the range of
+* the bivariate fields of the input grid.
+*
+* @section Case Tables
+* @subsection Marching tetrahedra with grey cases
+* In this filter, we have added one vertex classification in Marching Tetrahedra. The
+* reason is because we need a grey classification to ensure that surfaces coincident with
+* the boundary of the tetrahedra will also be included in the output. Given an iso-value,
+* each vertex on the tetrahedron can be classified into three types, including
+* equal-(G)rey, below-(W)hite or above-(B)lack the isovalue.
+* The notation of the classifications are represented as:
+* W or 0 --- below an iso-value
+* G or 1 --- equal an iso-value
+* B or 2 --- above an iso-value
+* The following illustration summarises all of the surface cases (Asterisk * is used to
+* highlight the outline of the triangle):
+* Case A (no triangles): 0000
+* No triangle is generated.
+*
+* W
+* ...
+* . . .
+* . . .
+* . .W. .
+* . . . .
+* W...........W
+*
+* Case B (one grey vertex): 0001, 0010, 0100, 1000
+* Only a vertex is kept, no triangle is generated.
+* W
+* ...
+* . . .
+* . . .
+* . .G. .
+* . . . .
+* W...........W
+*
+* Case C (one grey edge): 0011, 0101, 0110, 1001, 1010, 1100
+* Only an edge is kept, no triangle is generated.
+* G
+* ...
+* . . .
+* . . .
+* . .G. .
+* . . . .
+* W...........W
+*
+* Case D (standard triangle case): 0002, 0020, 0200, 2000
+* One triangle is generated
+* W W
+* ... ...
+* . . . . * .
+* . . . . *.* .
+* . .B. . -> . * B * .
+* . . . . . ** * ** .
+* W...........W W...........W
+*
+* Case E (one grey face): 0111, 1011, 1101, 1110
+* The triangle on one face of the tetra is generated.
+* G G
+* ... .**
+* . . . . * *
+* . . . -> . * *
+* . .G. . . .G* *
+* . . . . . . * *
+* W...........G W...........G
+*
+* Case F (triangle through vertex): 0012 0021 0102 0120 0201 0210
+* 1002 1020 1200 2001 2010 2100
+* A triangle connecting one of the tetra vertex is generated.
+* G G
+* ... .*.
+* . . . .*.*.
+* . . . -> . *.* .
+* . .B. . . *.B.* .
+* . . . . . * * * * .
+* W...........W W...........W
+*
+* Case G (grey tet - treat as empty): 1111
+* No triangle is generated.
+* G
+* ...
+* . . .
+* . . .
+* . .G. .
+* . . . .
+* G...........G
+*
+* Case H (triangle through edge): 0112 0121 0211 1012 1021 1102
+* 1120 1201 1210 2011 2101 2110
+* A triangle containing an edge of the tetra is generated.
+*
+* G G
+* ... ..*
+* . . . . . *
+* . . . . *. *
+* . . . . . *
+* . . . -> . * . *
+* . . W . . . . W . *
+* . . . . . * . *
+* . . . . . . * . *
+* B.............. G B...............G
+*
+* Case I (standard quad case): 0022 0202 0220 2002 2020 2200
+* A quand is generated, which can be split to two triangles.
+*
+* W W
+* ... ...
+* . . . . . .
+* . . . . . .
+* . . . * *. * *
+* . . . -> .* . *.
+* . . W . . . * . W . * .
+* . . . . . * * * * .
+* . . . . . . . .
+* B.............. B B..................B
+*
+* Case J (complement of Case E): 1112 1121 1211 2111
+* Case K (complement of Case F): 0122 0212 0221 1022 1202
+* 1220 2012 2021 2102 2120 2201 2210
+* Case L (complement of Case C): 1122 1212 1221 2112 2121 2211
+* Case M (complement of Case D): 0222 2022 2202 2220
+* Case N (complement of Case B): 1222 2122 2212 2221
+* Case O (complement of Case A): 2222
+*
+* @subsection Clipping cases of the base fiber surface
+* After generating the base fiber surface in each cell, we need a further clipping
+* process to obtain the accurate fiber surface. Clipping is based on classifying the
+* vertices of each triangle to several cases, which will be described in this section.
+* In order to keep things consistent, we assume that the vertices are ordered CCW, and
+* that edges are numbered according to the opposing vertex:
+*
+* v0
+* / \
+* e2 e1
+* / \
+* v1---e0---v2
+*
+* =======================================================================
+* There are six cases for clipping the base fiber surface (subject to the usual
+* symmetries & complementarity)
+*------------------------------------------------------------------------
+* Case A (No triangles): Cases 000 & 222
+* Here, the entire triangle is discarded
+*
+* 000(A): 0
+* . .
+* . .
+* . .
+* . .
+* . .
+* 0...........0
+*
+*------------------------------------------------------------------------
+* Case B (Point-triangle): Cases 001, 010, 100, 122, 212, 221
+* One vertex is kept, and a single triangle is extracted
+*
+* 001(B): 1
+* / \
+* / \
+* E-----E
+* . .
+* . .
+* 0...........0
+*
+*------------------------------------------------------------------------
+* Case C (Edge Quad): Cases 011, 101, 110, 112, 121, 211
+* Two vertices are kept, and a quad is extracted based on the edge
+*
+* 011(C): 1
+* /|\
+* / | \
+* / | E
+* / | / .
+* / |/ .
+* 1-----E.....0
+*
+*------------------------------------------------------------------------
+* Case D (Stripe): Cases 002, 020, 022, 200, 202, 220
+* No vertices are kept, but a stripe across the middle is generated
+*
+* 022(D): 2
+* . .
+* . E
+* . /|\
+* . / | E
+* . / |/ .
+* 2...E---E...0
+*
+*------------------------------------------------------------------------
+* Case E (Point-Stripe): Cases 012, 021, 102, 120, 201, 210
+* One vertex is kept, with a stripe through the triangle
+*
+* 021(E): 1
+* / \
+* E---E
+* .|\ |.
+* . | \ | .
+* . | \| .
+* 2...E---E...0
+*
+*------------------------------------------------------------------------
+* Case F (Entire): Case 111
+* All three vertices are kept, along with the triangle
+*
+* 111(F): 1
+* / \
+* / \
+* / \
+* / \
+* / \
+* 1-----------1
+*
+*
+* @section How to use this filter
+* Suppose we have a tetrahedral mesh stored in a vtkUnstructuredGrid, we call this
+* data set "inputData". This data set has two scalar arrays whose names are "f1"
+* and "f2" respectively. Assume the "inputPoly" is a valid input polygon. Given
+* these input above, this filter can be used as follows in c++ sample code:
+*
+* vtkSmartPointer<vtkFiberSurface> fiberSurface =
+* vtkSmartPointer<vtkFiberSurface>::New();
+* fiberSurface->SetInputData(0,inputData);
+* fiberSurface->SetInputData(1,inputPoly);
+* fiberSurface->SetField1("f1");
+* fiberSurface->SetField2("f2");
+* fiberSurface->Update();
+*
+* The filter output which can be called by "fiberSurface->GetOutput()", will be a
+* vtkPolyData containing the output fiber surfaces.
+*/
+
+#ifndef vtkFiberSurface_h
+#define vtkFiberSurface_h
+
+#include "vtkFiltersTopologyModule.h" // For export macro
+#include "vtkPolyDataAlgorithm.h"
+
+class VTKFILTERSTOPOLOGY_EXPORT vtkFiberSurface : public vtkPolyDataAlgorithm
+{
+public:
+ static vtkFiberSurface* New();
+ vtkTypeMacro(vtkFiberSurface, vtkPolyDataAlgorithm);
+ void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE;
+
+ /**
+ * Specify the first field name to be used in this filter.
+ */
+ void SetField1(char* fieldName);
+
+ /**
+ * Specify the second field name to be used in the filter.
+ */
+ void SetField2(char* fieldName);
+
+ /**
+ * This structure lists the vertices to use for the marching tetrahedra,
+ * Some of these vertices need to be interpolated, but others are the vertices of
+ * the tetrahedron already, and for this we need to store the type of vertex.
+ * bv_not_used: Symbol for an unused entry
+ * bv_vertex_*: Vertices on a tetrahedron
+ * bv_edge_*: Vertices on the edges of a tetrahedron
+ */
+ enum BaseVertexType
+ {
+ bv_not_used,
+ bv_vertex_0,
+ bv_vertex_1,
+ bv_vertex_2,
+ bv_vertex_3,
+ bv_edge_01,
+ bv_edge_02,
+ bv_edge_03,
+ bv_edge_12,
+ bv_edge_13,
+ bv_edge_23
+ };
+
+ /**
+ * After generating the base fiber surface in each cell, we need a further clipping
+ * process to obtain the accurate fiber surface. Clipping is based on classifying the
+ * vertices of each triangle, this structure lists the type of vertices to be used for
+ * the clipping triangles.
+ * Some of these vertices need to be interpolated, but others are the vertices of
+ * the triangle already, and for this we need to store the type of vertex.
+ * Moreover, vertices along the edges of the triangle may be interpolated either at
+ * parameter value 0 or at parameter value 1. In other words, there are three classes
+ * of vertex with three choices of each, with a total of nine vertices. We therefore
+ * store an ID code for each possibility as follows:
+ */
+ enum ClipVertexType
+ {
+ not_used,
+ vertex_0,
+ vertex_1,
+ vertex_2,
+ edge_0_parm_0,
+ edge_1_parm_0,
+ edge_2_parm_0,
+ edge_0_parm_1,
+ edge_1_parm_1,
+ edge_2_parm_1,
+ };
+
+protected:
+ vtkFiberSurface();
+ virtual int FillInputPortInformation(int port, vtkInformation* info) VTK_OVERRIDE;
+ virtual int RequestData(vtkInformation*, vtkInformationVector**, vtkInformationVector*) VTK_OVERRIDE;
+
+ // name of the input array names.
+ char* Fields[2];
+
+private:
+ vtkFiberSurface(const vtkFiberSurface&) VTK_DELETE_FUNCTION;
+ void operator=(const vtkFiberSurface&) VTK_DELETE_FUNCTION;
+};
+#endif
diff --git a/Infovis/Core/CMakeLists.txt b/Infovis/Core/CMakeLists.txt
index ece8bea78804e9d30e13cb05cc74783c7a1679a2..026fc9062aec6c4e1f4979e90a51699d675850ae 100644
--- a/Infovis/Core/CMakeLists.txt
+++ b/Infovis/Core/CMakeLists.txt
@@ -5,6 +5,7 @@ set(Module_SRCS
vtkArrayToTable.cxx
vtkCollapseGraph.cxx
vtkCollapseVerticesByArray.cxx
+ vtkContinuousScatterplot.cxx
vtkDataObjectToTable.cxx
vtkDotProductSimilarity.cxx
vtkExtractSelectedTree.cxx
diff --git a/Infovis/Core/Testing/Cxx/CMakeLists.txt b/Infovis/Core/Testing/Cxx/CMakeLists.txt
index e8488525317e903829eb7088ba6f9aebe724006a..14b9e40772446d2d491207d16cbcda48c7924c95 100644
--- a/Infovis/Core/Testing/Cxx/CMakeLists.txt
+++ b/Infovis/Core/Testing/Cxx/CMakeLists.txt
@@ -12,6 +12,7 @@ vtk_add_test_cxx(${vtk-module}CxxTests tests
ArrayTransposeMatrix.cxx,NO_VALID
TestArrayNorm.cxx,NO_VALID,NO_DATA
TestCollapseVerticesByArray.cxx,NO_VALID
+ TestContinuousScatterPlot.cxx,NO_VALID
TestDataObjectToTable.cxx,NO_VALID
TestExtractSelectedTree.cxx,NO_VALID
TestExtractSelectedGraph.cxx,NO_VALID
diff --git a/Infovis/Core/Testing/Cxx/TestContinuousScatterPlot.cxx b/Infovis/Core/Testing/Cxx/TestContinuousScatterPlot.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..d63ed183ff5b8de1c4ebd51a5c4d261dbc67c618
--- /dev/null
+++ b/Infovis/Core/Testing/Cxx/TestContinuousScatterPlot.cxx
@@ -0,0 +1,264 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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.
+
+=========================================================================*/
+
+// PURPOSE: test new vtkContinuousScatterplot filter
+#include <cstring>
+#include <sstream>
+#include <string>
+
+#include "vtkContinuousScatterplot.h"
+#include "vtkTestUtilities.h"
+#include <vtkActor.h>
+#include <vtkImageActor.h>
+#include <vtkImageData.h>
+#include <vtkImageMapper3D.h>
+#include <vtkImageWriter.h>
+#include <vtkJPEGWriter.h>
+#include <vtkPolyData.h>
+#include <vtkPolyDataMapper.h>
+#include <vtkRenderWindow.h>
+#include <vtkRenderWindowInteractor.h>
+#include <vtkRenderer.h>
+#include <vtkSmartPointer.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkXMLImageDataWriter.h>
+#include <vtkXMLUnstructuredGridReader.h>
+
+int TestContinuousScatterPlot(int argc, char* argv[])
+{
+ std::string outputString;
+ bool pass = true;
+
+ char* inputFile =
+ vtkTestUtilities::ExpandDataFileName(argc, argv, "Data/cube.vtu");
+
+ int ResX[5] = { 10, 20, 30, 40, 50 };
+ int ResY[5] = { 10, 20, 30, 40, 50 };
+
+ char fieldOneName[] = "f1";
+ char fieldTwoName[] = "f2";
+
+ int cmpIndex = 0;
+
+ /********************* Desired output arrays *********************/
+ std::string dataToCompare[5] = { "2,1,0,0,0,0,0,0,0,0,2,17,9,0,0,0,0,0,0,0,0,12,44,28,0,0,0,0,0,"
+ "0,0,0,40,81,64,0,0,0,0,0,0,0,0,42,137,112,0,0,0,0,0,0,0,0,48,"
+ "181,159,0,0,0,0,0,0,0,0,46,255,215,0,0,0,0,0,0,0,0,34,208,152,"
+ "0,0,0,0,0,0,0,0,18,157,55,0,0,0,0,0,0,0,0,5,61",
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,8,4,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,17,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,29,18,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,19,46,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,41,67,51,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,37,102,64,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,113,64,3,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,47,105,86,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,50,126,62,1,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,50,120,167,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,48,255,228,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,42,250,222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,35,225,197,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,26,183,157,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,18,134,108,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,11,97,56,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,47,14,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,1,16",
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,7,"
+ "3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,13,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,11,21,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "16,30,19,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,21,46,29,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,67,40,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,41,94,45,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,58,92,44,4,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,48,99,33,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,52,85,42,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,56,95,44,23,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,59,104,65,28,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,60,113,81,36,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,61,121,97,46,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,60,129,118,82,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,58,135,156,86,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,54,140,161,75,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,50,143,"
+ "150,44,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,146,237,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,255,225,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,29,209,191,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,22,168,153,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,126,113,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,11,90,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,54,"
+ "31,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,27,17,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0",
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,7,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,12,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,10,19,11,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,15,30,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,19,43,23,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,26,"
+ "55,31,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,31,64,29,2,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,39,74,29,4,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,49,78,36,7,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,55,76,46,12,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,42,86,36,28,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,44,70,53,36,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,47,60,64,34,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,1,28,64,63,37,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,1,27,68,71,48,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,27,116,86,60,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,74,123,120,78,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "75,144,143,34,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,"
+ "147,141,39,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,74,165,"
+ "156,163,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,71,170,252,"
+ "189,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,68,143,255,191,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,63,145,255,181,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,60,144,244,159,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,53,143,222,133,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,46,142,193,92,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,40,141,163,48,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,33,140,220,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,27,217,200,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,21,171,162,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,132,131,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,12,94,85,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,8,63,48,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,4,38,22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,"
+ "24,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,0",
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,2,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,5,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,8,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,15,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,11,22,11,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,15,29,15,1,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,19,37,"
+ "15,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,22,42,19,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,26,47,21,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,30,56,21,6,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,39,56,26,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,6,43,64,33,12,1,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,48,60,41,26,1,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,38,67,55,32,2,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,40,"
+ "71,48,38,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,5,41,58,57,37,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,6,42,62,69,53,5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,43,79,81,28,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,5,45,67,73,18,4,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,27,69,73,22,6,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,3,25,70,89,28,7,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,2,24,81,96,34,"
+ "14,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "2,22,84,103,46,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,1,21,85,88,55,22,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,21,66,95,65,27,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,20,111,102,88,44,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,65,111,109,94,43,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,62,111,117,95,35,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "59,111,124,89,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,55,110,132,193,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,53,109,255,196,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,48,106,231,178,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,104,214,158,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,38,102,191,"
+ "130,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,33,99,168,100,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,29,99,178,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,24,132,121,35,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,20,94,150,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,16,141,159,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,124,109,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,10,79,"
+ "77,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,7,55,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,5,37,28,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,3,22,13,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,11,7,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,"
+ "0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0" };
+
+ /************** output to error log ****************/
+ /****************************************************/
+
+ /****************************************************/
+
+ for (int i = 0; i < 5; i++)
+ {
+
+ // read .vtu file from the the command line
+ vtkSmartPointer<vtkXMLUnstructuredGridReader> mesh_reader =
+ vtkSmartPointer<vtkXMLUnstructuredGridReader>::New();
+ mesh_reader->SetFileName(inputFile);
+ mesh_reader->Update();
+
+ // Define the CSP filter.
+ vtkSmartPointer<vtkContinuousScatterplot> csp =
+ vtkSmartPointer<vtkContinuousScatterplot>::New();
+ csp->SetInputData(mesh_reader->GetOutput());
+ csp->SetField1(fieldOneName, ResX[cmpIndex]);
+ csp->SetField2(fieldTwoName, ResY[cmpIndex]);
+ csp->Update();
+
+ /************** Console output ****************/
+ {
+ std::stringstream ss;
+ int indexLocal = 0;
+
+ // Retrieve the entries from the image data
+ for (int z = 0; z < 1; z++)
+ {
+ for (int y = 0; y < ResY[cmpIndex]; y++)
+ {
+ for (int x = 0; x < ResX[cmpIndex]; x++)
+ {
+ double* pixel = static_cast<double*>(csp->GetOutput()->GetScalarPointer(x, y, z));
+ if (indexLocal > 0)
+ {
+ outputString.append(",");
+ }
+ ss << (int)(pixel[0] + 0.5);
+ outputString.append(ss.str());
+ ss.str(std::string());
+ ss.clear();
+ indexLocal++;
+ }
+ }
+ }
+
+ } // end of console output
+
+ /*********** comparing the strings ************/
+ if (strcmp(outputString.c_str(), dataToCompare[cmpIndex].c_str()) != 0)
+ {
+ pass = false;
+ { // writing to file
+ printf("\n\n/**************************************/");
+ printf("\n/********Test Unsuccessful*************/");
+ printf("\nInput Data: %s", (std::string(inputFile).substr(6, 49)).c_str());
+ printf("\nCS Plot test case %d x %d \n", ResX[cmpIndex], ResY[cmpIndex]);
+ printf("\nString to compare: %s", dataToCompare[cmpIndex].c_str());
+ printf("\nOutput String : %s", outputString.c_str());
+ printf("\n/**************************************/");
+ } // end of writing
+ }
+
+ cmpIndex++;
+ outputString = "";
+
+ } // end for
+
+ delete [] inputFile;
+
+ if (pass)
+ {
+ cout << "\nTest Successful!!!" << endl;
+ return EXIT_SUCCESS;
+ }
+
+ cout << "\nTest Unsuccessful." << endl;
+ return EXIT_FAILURE;
+}
diff --git a/Infovis/Core/vtkContinuousScatterplot.cxx b/Infovis/Core/vtkContinuousScatterplot.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5534b1baf45d6ecc3b5d173ccb4de9ec297378ad
--- /dev/null
+++ b/Infovis/Core/vtkContinuousScatterplot.cxx
@@ -0,0 +1,957 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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 "vtkContinuousScatterplot.h"
+
+#include "vtkCell.h"
+#include "vtkCellArray.h"
+#include "vtkCellData.h"
+#include "vtkCellType.h"
+#include "vtkDataSetAttributes.h"
+#include "vtkExecutive.h"
+#include "vtkFloatArray.h"
+#include "vtkIdList.h"
+#include "vtkIdTypeArray.h"
+#include "vtkImageData.h"
+#include "vtkInformation.h"
+#include "vtkInformationVector.h"
+#include "vtkMassProperties.h"
+#include "vtkMath.h"
+#include "vtkMergePoints.h"
+#include "vtkObjectFactory.h"
+#include "vtkPointData.h"
+#include "vtkPoints.h"
+#include "vtkPolyData.h"
+#include "vtkPolygon.h"
+#include "vtkSmartPointer.h"
+#include "vtkTriangleFilter.h"
+#include "vtkUnstructuredGrid.h"
+
+vtkStandardNewMacro(vtkContinuousScatterplot);
+
+// Data structure to store the fragment faces.
+// Each face of the fragment can be represented using a vtkIdList.
+typedef std::vector<vtkSmartPointer<vtkIdList> >* Polytope;
+
+//----------------------------------------------------------------------------
+vtkContinuousScatterplot::vtkContinuousScatterplot()
+{
+ // value for floating comparison. Suppose two floating values a and b:
+ // if fabs(a-b) <= dEpsilon : then a == b.
+ // if (a-b) > dEpsilon : then a > b.
+ // if (b-a) > dEpsilon : then b > a.
+ this->Epsilon = 1.0e-6;
+ // The number of output port is one.
+ this->SetNumberOfOutputPorts(1);
+ // Resolution of the output image is set to 100 as default.
+ this->ResX = this->ResY = 100;
+ // Create the vtkImageData object and pass to the output port.
+ vtkImageData* output = vtkImageData::New();
+ this->GetExecutive()->SetOutputData(0, output);
+ output->Delete();
+ this->Fields[0] = this->Fields[1] = NULL;
+}
+//----------------------------------------------------------------------------
+void vtkContinuousScatterplot::SetField1(char* nm, vtkIdType xRes)
+{
+ this->Fields[0] = nm;
+ this->ResX = xRes;
+}
+//----------------------------------------------------------------------------
+void vtkContinuousScatterplot::SetField2(char* nm, vtkIdType yRes)
+{
+ this->Fields[1] = nm;
+ this->ResY = yRes;
+}
+//----------------------------------------------------------------------------
+void vtkContinuousScatterplot::PrintSelf(ostream& os, vtkIndent indent)
+{
+ this->Superclass::PrintSelf(os, indent);
+}
+//----------------------------------------------------------------------------
+int vtkContinuousScatterplot::FillInputPortInformation(int, vtkInformation* info)
+{
+ info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkUnstructuredGrid");
+ return 1;
+}
+//----------------------------------------------------------------------------
+int vtkContinuousScatterplot::FillOutputPortInformation(int, vtkInformation* info)
+{
+ info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkImageData");
+ return 1;
+}
+//----------------------------------------------------------------------------
+int vtkContinuousScatterplot::RequestData(
+ vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
+{
+ // get the input and output ports information.
+ vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
+ vtkInformation* outInfo = outputVector->GetInformationObject(0);
+
+ // get the input data set, which is required to be a vtkUnstructuredGrid.
+ vtkUnstructuredGrid* input =
+ vtkUnstructuredGrid::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
+ // get the output data set, which is required to be a vtkImageData.
+ vtkImageData* output = vtkImageData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
+
+ // get the number of points and cells of the input grid.
+ vtkIdType numCells = input->GetNumberOfCells();
+ vtkIdType numPts = input->GetNumberOfPoints();
+
+ // get the dataset statistics and check if it is not an empty grid.
+ if (numCells < 1 || numPts < 1)
+ {
+ vtkErrorMacro("No input data.");
+ return 1;
+ }
+
+ // check if the output image resolution is a valid number.
+ if (this->ResX <= 0 || this->ResY <= 0)
+ {
+ vtkErrorMacro("The resolution of the output image has to be a positive number.");
+ return 1;
+ }
+
+ // check if the names of the input scalar fields are specified.
+ if (!this->Fields[0] || !this->Fields[1])
+ {
+ vtkErrorMacro("At least two fields need to be specified.");
+ return 1;
+ }
+
+ // Input point data, which should include the arrays defining the range space.
+ vtkDataSetAttributes* inPD = input->GetPointData();
+ // current data array of the given array names from the user input.
+ vtkSmartPointer<vtkDataArray> array;
+ // scalar range in two fields, where fieldInterval[0] = f1_max - f1_min,
+ // fieldInterval[1] = f2_max - f2_min
+ float fieldInterval[2];
+
+ // fragment width of two fields, where fragWidth[0] = fieldInterval[0] / resX,
+ // fragWidth[0] = fieldInterval[1] / resY, (resX, resY) is the output image
+ // resolution.
+ float fragWidth[2];
+
+ // Collect the first scalar field based on array names,
+ array = inPD->GetArray(this->Fields[0]);
+ // confirming that this field is present in the input.
+ if (array)
+ {
+ // Collect field ranges for later use.
+ fieldInterval[0] = array->GetRange()[1] - array->GetRange()[0];
+ // Interval between cutting planes in this field.
+ fragWidth[0] = (float)fieldInterval[0] / this->ResX;
+ }
+ else
+ {
+ vtkErrorMacro("Array not found in input point data: " << this->Fields[0]);
+ return 1;
+ }
+
+ // Collect the second field data array.
+ array = inPD->GetArray(this->Fields[1]);
+ // confirming that this field is present in the input.
+ if (array)
+ {
+ // The range interval of the field.
+ fieldInterval[1] = array->GetRange()[1] - array->GetRange()[0];
+ // Interval between cutting planes in this field.
+ fragWidth[1] = (float)fieldInterval[1] / this->ResY;
+ }
+ else
+ {
+ vtkErrorMacro("Array not found in input point data: " << this->Fields[1]);
+ return 1;
+ }
+
+ // Devide the tetrahedron into four faces. The index of each face.
+ const int tetTemplate[4][3] = { { 0, 1, 2 }, { 0, 1, 3 }, { 0, 2, 3 }, { 1, 2, 3 } };
+
+ // fragments of current cell: each cell is placed into the inputQ,
+ // then fresh fragments from current slice put into outputQ.
+ std::vector<Polytope> inputQ = std::vector<Polytope>();
+ std::vector<Polytope> outputQ = std::vector<Polytope>();
+
+ // structure for storing the vertices of cutting planes used to subdivide the cell.
+ vtkSmartPointer<vtkIdList> cut = vtkSmartPointer<vtkIdList>::New();
+
+ // Initially the points in the cutting plane are not ordered. These parameters are used
+ // to compute convex hull for points lying on the cutting plane. The sorting process
+ // utilizes the angular information between the plane edges.
+ float theta[100];
+ vtkIdType index[100];
+ // number of points in the current cutting plane.
+ int nrPoints;
+ float base[3], vec[3];
+ float baseNorm, ang;
+ vtkIdType pnt, pnt0, pnt1;
+
+ // min and max scalar value of the current field in a tetrahedron
+ float maxCell, minCell;
+ // min and max scalar values of the current field in the whole domain.
+ float minField = 0;
+ float maxField = 0;
+
+ // first and last cutting planes' field values. This threshold is used to traverse
+ // through all of the cutting planes in the range.
+ float initThreshold, lastThreshold;
+
+ // variables for walking around a cell face.
+ vtkIdType newPointId, thisPointId, prevPointId;
+ // number of the vertices in the current cell face.
+ int nrFaceIds;
+
+ // For each cutting plane, determine its intersecting position along the current edge.
+ int mask;
+
+ // scalar values of the end points of the current edge.
+ float thisScalar, prevScalar;
+
+ // if the cutting plane intersects in the middle of the given edge, these parameters
+ // are used to compute these intersecting points using linear interpolation
+ double delta, t, p[3], p0[3], p1[3], p2[3];
+
+ // fragment: list of faces forming the fragments in the current cell.
+ // residual: list of faces forming the rest of the cell (except fragments).
+ // working: collect the residual faces for the next iteration of subdivision.
+ Polytope fragment = NULL, residual = NULL, working = NULL;
+
+ // structure for storing the current framgent vertices in each cell face.
+ vtkSmartPointer<vtkIdList> fragmentFace = NULL;
+ // structure for storing the vertices which are not belonging to the current
+ // fragment in each cell face.
+ vtkSmartPointer<vtkIdList> residualFace = NULL;
+
+ // reading/writing cells from input, output and cell arrays.
+ vtkSmartPointer<vtkIdList> cell = NULL;
+ // reading/writing faces from input, output and cell arrays.
+ vtkSmartPointer<vtkIdList> face = NULL;
+
+ // Each fragment can be mapped to a 2-D bin based on its range values. The density
+ // value in each bin equals to the total geometric volume of the fragments in this bin.
+ // The following structure creates and initialise such a 2-D bin with a resolution
+ // ResX * ResY.
+ // float imageBin[this->ResX][this->ResY];
+ float** imageBin = new float*[this->ResX];
+ for (int xIndex = 0; xIndex < this->ResX; ++xIndex)
+ {
+ imageBin[xIndex] = new float[this->ResY];
+ for (int yIndex = 0; yIndex < this->ResY; ++yIndex)
+ {
+ imageBin[xIndex][yIndex] = 0.0;
+ }
+ }
+ // maximal volume value in the output image bin.
+ float maxBinSize = 0;
+
+ // The VTK filter for computing the fragment geometric volume. This filter only
+ // accepts triangular mesh as input, while our fragments are stored as a polygonal
+ // mesh. Therefore, we need a mesh conversion.
+ vtkSmartPointer<vtkMassProperties> volume = vtkSmartPointer<vtkMassProperties>::New();
+ // This VTK filter convert a polygonal mesh to a triangular mesh.
+ vtkSmartPointer<vtkTriangleFilter> triangleFilter = vtkSmartPointer<vtkTriangleFilter>::New();
+ // geometric volume of the current polyhedral fragments.
+ float fragVolume;
+
+ // mapped index of a fragment into the 2-D bin.
+ vtkIdType binIndexFirst = 0, binIndexSecond = 0;
+
+ // id for the newly interpolated fragment vertices. Since we only need to update point
+ // structure, therefore this id is not of interest and can be ignored for the moment.
+ vtkIdType ignored = 0;
+ // point coordinates in the input grid
+ vtkPoints* inputGridPoints = input->GetPoints();
+ // point coordinates in the output fragments
+ vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
+
+ // Each fragment is composed of a number of points. This data structure contains
+ // scalar values for all these points.
+ vtkSmartPointer<vtkDataSetAttributes> newPointsPD = vtkSmartPointer<vtkDataSetAttributes>::New();
+
+ // locator for inserting new interpolated points into the structure. This search
+ // structure ensures that redundant points will not be added to the point collection.
+ vtkSmartPointer<vtkMergePoints> pointLocator = vtkSmartPointer<vtkMergePoints>::New();
+
+ // Each fragment will have a range value. This array records the range values
+ // for all of the fragments. Since we iteratively subdivide the cell based on
+ // two fields, therefore the first n elements in this array record the
+ // range values of the fragments in the first field subdivision, and the following m
+ // elements record the range values in the second field subdivision.
+ vtkSmartPointer<vtkFloatArray> fragScalar = vtkSmartPointer<vtkFloatArray>::New();
+
+ // Each fragment will have range values from two fields.
+ float fragRangeValue[2];
+ // number of fragments in the first field subdivision.
+ size_t firstFragNr = 0;
+
+ // data structure containing scalar values of the vertices of current tetrahedron of
+ // the input grid.
+ vtkSmartPointer<vtkDataSetAttributes> tetraPD = vtkSmartPointer<vtkDataSetAttributes>::New();
+
+ // scalar values of two fields in the current tetrahedron of the input grid.
+ // these two arrays are added to the tetrahedron point data.
+ vtkSmartPointer<vtkFloatArray> tetraF1 = vtkSmartPointer<vtkFloatArray>::New(),
+ tetraF2 = vtkSmartPointer<vtkFloatArray>::New();
+ tetraF1->SetNumberOfComponents(1);
+ tetraF1->SetNumberOfTuples(4);
+ tetraF2->SetNumberOfComponents(1);
+ tetraF2->SetNumberOfTuples(4);
+ // add the two scalar arrays to the point data of the tetrahedron
+ tetraPD->AddArray(tetraF1);
+ tetraPD->AddArray(tetraF2);
+
+ // structure containing the polygonal faces of polyhedral fragments
+ vtkSmartPointer<vtkPolyData> polyhedra = vtkSmartPointer<vtkPolyData>::New();
+
+ // estimate the total number of cells of the fragments in each input cell.
+ // consider edge can be maximally subdivided to resX number
+ // of fragments for the first field and resY number of fragments for the second field.
+ // In total, there will be maximal resX * resY number of new points in each edge.
+ int estOutputPointSize = this->ResX * this->ResY * 4;
+ // Allocate the memory for the framgent points.
+ newPoints->Allocate(estOutputPointSize);
+
+ // main loop ...
+ // For each tetrahedron in a gird
+ for (vtkIdType tetraIndex = 0; tetraIndex < input->GetNumberOfCells(); tetraIndex++)
+ {
+ // current tetrahedron vertex list.
+ cell = vtkSmartPointer<vtkIdList>::New();
+ input->GetCellPoints(tetraIndex, cell);
+
+ // Test if the current cell is a tetrahedron or not. If not, ignore this cell.
+ if (input->GetCellType(tetraIndex) != VTK_TETRA || cell->GetNumberOfIds() != 4)
+ {
+ vtkWarningMacro("Current cell " << tetraIndex << " is not of a tetrahedron type.");
+ continue;
+ }
+
+ // initialise the point structure for storing fragment vertices.
+ newPoints->Reset();
+ // initialise the search structure for new points insertion.
+ pointLocator->InitPointInsertion(newPoints, input->GetCell(tetraIndex)->GetBounds());
+
+ // initialise data structure containing the scalar values for the fragment vertices.
+ newPointsPD->Initialize();
+ // Set the storage for interpolating the field values of the fragment vertices.
+ newPointsPD->InterpolateAllocate(tetraPD, estOutputPointSize, ignored, false);
+ newPointsPD->CopyScalarsOn();
+
+ // initialise data structure containing the scalar values of the whole fragment.
+ fragScalar->Initialize();
+ fragScalar->Allocate(this->ResX * this->ResY);
+ // two components are needed to store the bivariate fields of the framgent.
+ fragScalar->SetNumberOfComponents(2);
+
+ // Initialise the scalar values in this tetrahedral cell.
+ for (vtkIdType cellIndex = 0; cellIndex < cell->GetNumberOfIds(); cellIndex++)
+ {
+ int pointId = cell->GetId(cellIndex);
+ pointLocator->InsertNextPoint(inputGridPoints->GetPoint(pointId));
+ tetraF1->SetComponent(
+ cellIndex, 0, inPD->GetArray(this->Fields[0])->GetComponent(pointId, 0));
+ tetraF2->SetComponent(
+ cellIndex, 0, inPD->GetArray(this->Fields[1])->GetComponent(pointId, 0));
+ }
+
+ // The scalar values of the framgent points are based on the interpolation of the
+ // point data of the tetrahedral cell (tetraPD).
+ for (vtkIdType cellPDIndex = 0; cellPDIndex < tetraPD->GetNumberOfTuples(); cellPDIndex++)
+ {
+ newPointsPD->CopyData(tetraPD, cellPDIndex, cellPDIndex);
+ }
+
+ /*
+ cout << "ARRAY SUMMARY:\n" << endl;
+ for (int i = 0; i < tetraPD->GetNumberOfTuples(); i++)
+ {
+ cout << i << "\t" << inPD->GetArray(this->Fields[0])->GetComponent(i,0)
+ << "\t" << inPD->GetArray(this->Fields[1])->GetComponent(i,0)
+ << "\t" << newPointsPD->GetArray(0)->GetComponent(i,0)
+ << "\t" << newPointsPD->GetArray(1)->GetComponent(i,0) << endl;
+
+ }
+
+ */
+ // Get the next cell from input, and place into the working queue.
+ // We place into outputQ as each field takes the output of the
+ // last step as its input, swapping queues BEFORE processing.
+ // OutputQ : a list of faces of the output framgent
+ Polytope ptp = new std::vector<vtkSmartPointer<vtkIdList> >();
+ outputQ.push_back(ptp);
+
+ // min and max range value of the current cell
+ minCell = maxField;
+ maxCell = minField;
+
+ // for each point in the cell
+ for (int fnr = 0; fnr < 4; fnr++)
+ {
+ // get the faces of the cell and push into outputQ structure.
+ face = vtkSmartPointer<vtkIdList>::New();
+ face->SetNumberOfIds(3);
+ for (int pnr = 0; pnr < 3; pnr++)
+ {
+ face->SetId(pnr, tetTemplate[fnr][pnr]);
+ }
+ outputQ[0]->push_back(face);
+ }
+
+ // For each scalar field:
+ for (size_t fieldNr = 0; fieldNr < 2; fieldNr++)
+ {
+ // Swap role of outputQ and inputQ
+ std::swap(outputQ, inputQ);
+ outputQ.clear();
+
+ // If this is the second round subdivision, then we need to record the number
+ // of fragments produced in the first round. This number is used to locate
+ // the range values of the output fragment.
+ if (fieldNr == 1)
+ {
+ firstFragNr = inputQ.size();
+ }
+
+ // min value of the whole field
+ minField = inPD->GetArray(this->Fields[fieldNr])->GetRange()[0];
+ // max value of the whole field
+ maxField = inPD->GetArray(this->Fields[fieldNr])->GetRange()[1];
+ // Initialize the value
+ minCell = maxField;
+ maxCell = minField;
+
+ ////cout << "tet " << tetraIndex << ", field " << fieldNr << ", pnr: " <<
+ ///newPointsPD->GetNumberOfTuples() << endl;
+ ////cout << "min/max init: " << minCell << ", " << maxCell << endl;
+
+ // obtain the minimal and maximal scalar values of the cell.
+ for (int pnr = 0; pnr < newPointsPD->GetNumberOfTuples(); pnr++)
+ {
+ double fval = newPointsPD->GetArray((int)fieldNr)->GetComponent(pnr, 0);
+
+ ////cout << " fval: " << fval << endl;
+ if (maxCell < fval)
+ {
+ maxCell = fval;
+ }
+ if (minCell > fval)
+ {
+ minCell = fval;
+ }
+ }
+
+ ////cout << "D" << endl;
+ ////cout << "Cell min/max: " << minCell << " " << maxCell << endl;
+ ////cout << "field min/max[0] " << minField << " " << maxField << endl;
+ ////cout << "field widths " << fragWidth[0] << " " << fragWidth[1] << endl;
+
+ // in each field, the smallest threshold of cutting plane to start with.
+ // since each field is sliced uniformly, in other words, the interval between
+ // fragment in each field is a constant. Therefore in each cell, the smallest
+ // threshold to start with, should be the first fragment value of the field which
+ // is greater or equal to the minimal range value of the cell.
+ initThreshold =
+ minField + (1 + floor((minCell - minField) / fragWidth[fieldNr])) * fragWidth[fieldNr];
+
+ // Iterate through the faces of the current mesh.
+ // 1. The first field is initially traversed and processed on the input grid.
+ // The inputQ should be the initial tetrahedron mesh. And the outputQ should be
+ // the collection of computed fragments with respect to the first field.
+ // 2. The second field is then processed based on the outputQ of step 1. The inputQ
+ // now becomes the fragments in the first field subdivision. The OutputQ
+ // in this step should be the fragments of the bivariate fields.
+ // Fragment: array records the faces of the fragments.
+ // Residual: array records the faces of rest of the cell in addition to fragments.
+ // Cut: array records the vertices of the current cutting plane.
+ // For each faces of the input mesh.
+ for (size_t cp = 0; cp < inputQ.size(); cp++)
+ {
+ // working[0] : faces of current input mesh.
+ working = inputQ[cp];
+
+ // initialise the first and last cutting planes in this cell.
+ lastThreshold = initThreshold;
+
+ // Traverse from the minimal to the maximal scalar values in a cell, every time
+ // the threshold is increased by one fragmentWidth
+ for (double threshold = initThreshold; threshold < maxCell; threshold += fragWidth[fieldNr])
+ {
+ // Initialise framgent face structure for the current cutting plane.
+ if (fragment)
+ {
+ delete fragment;
+ }
+ if (residual)
+ {
+ delete residual;
+ }
+ fragment = new std::vector<vtkSmartPointer<vtkIdList> >();
+ residual = new std::vector<vtkSmartPointer<vtkIdList> >();
+
+ // Create the new cutting plane.
+ cut = vtkSmartPointer<vtkIdList>::New();
+
+ // Effectively, we start processing a new cell at this point.
+ for (std::vector<vtkSmartPointer<vtkIdList> >::iterator faceIt = working->begin();
+ faceIt != working->end(); faceIt++)
+ {
+ fragmentFace = vtkSmartPointer<vtkIdList>::New();
+ residualFace = vtkSmartPointer<vtkIdList>::New();
+
+ // number of points in the current cell face
+ nrFaceIds = (*faceIt)->GetNumberOfIds();
+
+ // get the previous point id in the face
+ prevPointId = (*faceIt)->GetId(nrFaceIds - 1);
+ // the scalar value of the previous point in the face
+ prevScalar = newPointsPD->GetArray(
+ (int)fieldNr)->GetComponent(prevPointId, 0);
+
+ // Walk around the edge, comparing the range values between the current
+ // cutting plane and the edge end points. Classify the each end point of the
+ // edge into three classes based on the comparison result. Three classes
+ // includes:
+ // fragmentFace: if current edge point belongs to the fragment.
+ // cut: if current edge point belongs to the cutting plane.
+ // residual: if current edge point belongs to neither of the above classes.
+ // For each edge
+ for (int i = 0; i < nrFaceIds; i++)
+ {
+ // get the current point Id in the face
+ thisPointId = (*faceIt)->GetId(i);
+ // get scalar value of the current point
+ thisScalar = newPointsPD->GetArray(
+ (int)fieldNr)->GetComponent(thisPointId, 0);
+
+ ////cout << ">>> " << thisPointId << " " << thisScalar << " " << prevPointId << " "
+ ///<< prevScalar << endl;
+
+ // zero bitweise or to any value equals that value
+ // threshold = minField + n * sw
+ mask = 0;
+ // thisScalar < threshold (0001)
+ if (threshold - thisScalar > this->Epsilon)
+ mask |= 1;
+ // thisScalar > threshold (0010)
+ if (thisScalar - threshold > this->Epsilon)
+ mask |= 2;
+ // prevScalar < threshold (0100)
+ if (threshold - prevScalar > this->Epsilon)
+ mask |= 4;
+ // prevScalar > threshold (1000)
+ if (prevScalar - threshold > this->Epsilon)
+ mask |= 8;
+ // thisScalar == threshold (0011)
+ if (fabs(thisScalar - threshold) <= this->Epsilon)
+ mask |= 3;
+ // prevScalar == threshold (1100)
+ if (fabs(prevScalar - threshold) <= this->Epsilon)
+ mask |= 12;
+ switch (mask)
+ {
+ // threshold == thisScalar
+ case 3:
+ case 7:
+ // threshold == thisScalar,
+ // threshold < preScalar
+ case 11:
+ // threshold == thisScalar
+ // threshold == preScalar
+ case 15:
+ // fragment face array insert current point
+ fragmentFace->InsertNextId(thisPointId);
+ // residual face array insert current point
+ residualFace->InsertNextId(thisPointId);
+ // if the current point is not in the cut array
+ if (cut->IsId(thisPointId) < 0)
+ {
+ cut->InsertNextId(thisPointId);
+ }
+ break;
+ // threshold == thisScalar
+ // threshold > preScalar
+ case 5:
+ // threshold > thisScalar
+ // threshold == prevScalar
+ case 13:
+ fragmentFace->InsertNextId(thisPointId);
+ break;
+ // threshold < thisScalar
+ // threshold < prevScalar
+ case 10:
+ // threshold < thisScalar
+ // threshold == prevScalar
+ case 14:
+ residualFace->InsertNextId(thisPointId);
+ break;
+ // threshold > thisScalar
+ // threshold < preScalar
+ case 9:
+ // Moving to this point crosses the threshold hi-lo
+ // Insert interpolated point into both.
+ delta = prevScalar - thisScalar;
+ t = (threshold - thisScalar) / delta;
+ // * PREV ------- T -------- THIS *
+ newPoints->GetPoint(thisPointId, p1);
+ newPoints->GetPoint(prevPointId, p2);
+ for (int j = 0; j < 3; j++)
+ {
+ p[j] = p1[j] + t * (p2[j] - p1[j]);
+ }
+ if (pointLocator->InsertUniquePoint(p, newPointId))
+ {
+ newPointsPD->InterpolateEdge(
+ newPointsPD, newPointId, thisPointId, prevPointId, t);
+ }
+ fragmentFace->InsertNextId(newPointId);
+ fragmentFace->InsertNextId(thisPointId);
+ residualFace->InsertNextId(newPointId);
+ // We have found a cut point, add it if first visit.
+ if (cut->IsId(newPointId) < 0)
+ {
+ cut->InsertNextId(newPointId);
+ }
+ break;
+ // threshold < thisScalar
+ // threshold > preScalar
+ case 6:
+ // this >, prev <
+ // Moving to this point crosses the threshold lo-hi
+ // Insert interpolated point into both.
+ delta = thisScalar - prevScalar;
+ t = (threshold - prevScalar) / delta;
+ newPoints->GetPoint(thisPointId, p1);
+ newPoints->GetPoint(prevPointId, p2);
+ for (int j = 0; j < 3; j++)
+ {
+ p[j] = p2[j] + t * (p1[j] - p2[j]);
+ }
+ if (pointLocator->InsertUniquePoint(p, newPointId))
+ {
+ newPointsPD->InterpolateEdge(
+ newPointsPD, newPointId, prevPointId, thisPointId, t);
+ }
+ fragmentFace->InsertNextId(newPointId);
+ residualFace->InsertNextId(newPointId);
+ residualFace->InsertNextId(thisPointId);
+ if (cut->IsId(newPointId) < 0)
+ {
+ cut->InsertNextId(newPointId);
+ }
+ break;
+ default:
+ vtkErrorMacro("Incomparable scalars " << prevScalar << ", " << thisScalar << ", "
+ << threshold);
+ } // switch mask
+ prevPointId = thisPointId;
+ prevScalar = thisScalar;
+ } // for-each edge
+ // Output fragment and residue into new cells, as appropriate.
+ // Test if the fragment and residual faces are well defined.
+ if (fragmentFace->GetNumberOfIds() > 2)
+ {
+ fragmentFace->Squeeze();
+ fragment->push_back(fragmentFace);
+ }
+ if (residualFace->GetNumberOfIds() > 2)
+ {
+ residualFace->Squeeze();
+ residual->push_back(residualFace);
+ }
+ } // for each face.
+
+ // We need to compute the face defined by the cut points.
+ // We cannot guarantee that points in the cut-list are ordered wrt
+ // polygon boundary, so we recompute an order by effectively working
+ // the convex hull.
+ // The cut-list is added to the framgent and residual array.
+ if (cut->GetNumberOfIds() > 2)
+ {
+ nrPoints = cut->GetNumberOfIds();
+ pnt0 = cut->GetId(0);
+ pnt1 = cut->GetId(1);
+ // 2. Compute vector from p[0] to p[1].
+ newPoints->GetPoint(pnt0, p0);
+ newPoints->GetPoint(pnt1, p1);
+ for (int i = 0; i < 3; i++)
+ {
+ base[i] = p1[i] - p0[i];
+ }
+ baseNorm = vtkMath::Norm(base);
+ theta[0] = 0.0;
+ index[0] = pnt1;
+ // 3. For the remaining points p, compute angle between p-p0 and base.
+ for (int i = 2; i < nrPoints; i++)
+ {
+ newPoints->GetPoint(cut->GetId(i), p1);
+ for (int j = 0; j < 3; j++)
+ {
+ vec[j] = p1[j] - p0[j];
+ }
+ theta[i - 1] = acos(vtkMath::Dot(base, vec) / (baseNorm * vtkMath::Norm(vec)));
+ index[i - 1] = cut->GetId(i);
+ }
+ // 4. Sort angles.
+ for (int j = 1; j < nrPoints - 1; j++)
+ {
+ ang = theta[j];
+ pnt = index[j];
+ int i = j - 1;
+ while (i >= 0 && theta[i] > ang)
+ {
+ theta[i + 1] = theta[i];
+ index[i + 1] = index[i];
+ i--;
+ }
+ theta[i + 1] = ang;
+ index[i + 1] = pnt;
+ }
+ cut->Reset();
+ cut->InsertNextId(pnt0);
+ for (int i = 0; i < nrPoints - 1; i++)
+ {
+ cut->InsertNextId(index[i]);
+ }
+ fragment->push_back(cut);
+ residual->push_back(cut);
+ } // Generate cut plane. // Generate cut plane.
+
+ // OUTPUT PHASE -----------------------------------------
+ // Have completed traversing each face.
+ // Now update output with new fragment.
+ // Iterate through each edge and record for each fragment, which are the
+ // divided points included.
+ // ------------------------------------------------------
+ if (fragment->size() > 3)
+ {
+ // add the current framgent to the outpuQ structure.
+ outputQ.push_back(fragment);
+ fragment = NULL;
+ // set threshold at which this fragment was created
+ if (fieldNr)
+ {
+ // if this is the second field subdivision, we need to retrieve
+ // the scalar value from the first field subdivision first. And insert the
+ // retrieved first field value together with the current second field value
+ // into the fragment scalar array.
+ float firstFieldRange = fragScalar->GetComponent((vtkIdType)cp, 0);
+ fragScalar->InsertNextTuple2(firstFieldRange, threshold);
+ }
+ else
+ {
+ // if this is the first field subdivisions, just insert the scalar value
+ // into the array.
+ fragScalar->InsertNextTuple2(threshold, 0);
+ }
+ }
+ else
+ {
+ // Remove any partial fragments.
+ while (fragment->size() > 0)
+ {
+ fragment->pop_back();
+ }
+ }
+ // Faces defining the next working polyhedron are in the residual array.
+ std::swap(working, residual);
+ // Clear out anything left in residual
+ while (residual->size() > 0)
+ {
+ residual->pop_back();
+ }
+ lastThreshold += fragWidth[fieldNr]; // Track final threshold
+ } // for each threshold
+
+ // If we are left with a residual, add it as a fragment.
+ // However, note that final step in loop is to swap residual
+ // with working, so need to look in the working list.
+ if (working->size() > 3)
+ {
+ outputQ.push_back(working);
+ // set threshold at which this fragment was created.
+ if (fieldNr)
+ {
+ // if this is the second field subdivision, we need to retrieve
+ // the scalar value from the first field subdivision first. And insert the
+ // retrieved first field value together with the current second field value
+ // into the fragment scalar array.
+ float firstFieldRange = fragScalar->GetComponent((vtkIdType)cp, 0);
+ fragScalar->InsertNextTuple2(firstFieldRange, lastThreshold);
+ }
+ else
+ {
+ // if this is the first field subdivisions, just insert the scalar value
+ // into the array.
+ fragScalar->InsertNextTuple2(lastThreshold, 0);
+ }
+ }
+ else
+ {
+ while (working->size() > 0)
+ {
+ working->pop_back();
+ }
+
+ if(working)
+ {
+ delete working;
+ working = nullptr;
+ }
+ }
+ } // for each fragment of cell: faces, edges.
+ // CutIds are the edges with current dividing fragments.
+ } // for each field
+
+ // OUTPUT PHASE: ----------------------------------------------
+ // Output generated fragments into main output dataset.
+ // ------------------------------------------------------------
+ // For each output framgent, we compute its geometric volume and aggregate over
+ // the cells.
+ for (size_t co = 0; co < outputQ.size(); co++)
+ {
+ // The current fragment needs to be converted to a polygonal mesh.
+ // Array for recording the vertices of the polygonal mesh.
+ polyhedra->Initialize();
+ polyhedra->Allocate((vtkIdType)outputQ[co]->size());
+
+ // for each face of the fragment
+ vtkSmartPointer<vtkIdList> poly = vtkSmartPointer<vtkIdList>::New();
+ for (std::vector<vtkSmartPointer<vtkIdList> >::iterator fc = outputQ[co]->begin();
+ fc != outputQ[co]->end(); fc++)
+ {
+ poly->Reset();
+ for (int pnr = 0; pnr < (*fc)->GetNumberOfIds(); pnr++)
+ {
+ poly->InsertNextId((*fc)->GetId(pnr));
+ }
+ // insert the polygon face.
+ polyhedra->InsertNextCell(VTK_POLYGON, poly);
+ }
+ polyhedra->SetPoints(newPoints);
+
+ // convert the polygon faces to triangular faces.
+ triangleFilter->SetInputData(polyhedra);
+ triangleFilter->Update();
+
+ if (triangleFilter->GetOutput()->GetNumberOfCells() > 0)
+ {
+ // Compute the volume of the fragment
+ volume->SetInputData(triangleFilter->GetOutput());
+ volume->Update();
+ fragVolume = volume->GetVolume();
+ }
+ else
+ {
+ fragVolume = 0;
+ }
+
+ // Range values for the current fragment.
+ fragRangeValue[0] = fragScalar->GetComponent((vtkIdType)(firstFragNr + co), 0);
+ fragRangeValue[1] = fragScalar->GetComponent((vtkIdType)(firstFragNr + co), 1);
+
+ // Map the current fragment into the 2-Dimensional bin based on its range values.
+ binIndexFirst = (int)(this->ResX - 1) *
+ (fragRangeValue[0] - inPD->GetArray(this->Fields[0])->GetRange()[0]) / fieldInterval[0];
+ binIndexSecond = (int)(this->ResY - 1) *
+ (fragRangeValue[1] - inPD->GetArray(this->Fields[1])->GetRange()[0]) / fieldInterval[1];
+
+ ////cout << "biF: " << binIndexFirst << "\tbiS: " << binIndexSecond << "\t" << fragVolume <<
+ ///endl;
+
+ // aggregate the fragment volumes in each bin
+ if (binIndexFirst >= 0 && binIndexFirst < this->ResX && binIndexSecond >= 0 &&
+ binIndexSecond < this->ResY)
+ {
+ imageBin[binIndexFirst][binIndexSecond] += fragVolume;
+ }
+
+ // finding the largest volume in a bin.
+ if (imageBin[binIndexFirst][binIndexSecond] > maxBinSize)
+ {
+ maxBinSize = imageBin[binIndexFirst][binIndexSecond];
+ }
+
+ // Clear faces from current polytope in output queue.
+ while (outputQ[co]->size() > 0)
+ {
+ outputQ[co]->pop_back();
+ }
+ } // end for loop in outputQ
+
+ // Release the memory for the local structure.
+ while (outputQ.size() > 0)
+ {
+ Polytope plocal;
+ plocal = outputQ.back();
+ outputQ.pop_back();
+ delete plocal;
+ }
+ outputQ.clear();
+ } // For each cell
+
+ if (fragment)
+ {
+ delete fragment;
+ fragment = nullptr;
+ }
+ if (residual)
+ {
+ delete residual;
+ residual = nullptr;
+ }
+
+
+ // Create the output image data.
+ output->SetExtent(0, this->ResX - 1, 0, this->ResY - 1, 0, 0);
+ output->SetOrigin(0, 0, 0);
+ output->SetSpacing(1, 1, 1);
+ output->AllocateScalars(VTK_DOUBLE, 1);
+
+ // A scalar array is attached onto the output vtkImageData.
+ // This array records the total volume of the fragments in each bin.
+ vtkSmartPointer<vtkFloatArray> volumeArray = vtkSmartPointer<vtkFloatArray>::New();
+ volumeArray->SetName("volume");
+ volumeArray->SetNumberOfComponents(1);
+ volumeArray->SetNumberOfTuples(this->ResX * this->ResY);
+
+ // Pixel density are computed by the fragment volume.
+ for (int xIndex = 0; xIndex < this->ResX; xIndex++)
+ {
+ for (int yIndex = 0; yIndex < this->ResY; yIndex++)
+ {
+ float size = imageBin[xIndex][yIndex];
+ int idx = xIndex * this->ResY + yIndex;
+ volumeArray->SetComponent(idx, 0, size);
+ double* pixel = static_cast<double*>(output->GetScalarPointer(xIndex, yIndex, 0));
+ if (imageBin[xIndex][yIndex] > 0)
+ {
+ pixel[0] = 255.0 * imageBin[xIndex][yIndex] / maxBinSize;
+ }
+ else
+ {
+ pixel[0] = 0;
+ }
+ }
+ }
+ output->GetPointData()->AddArray(volumeArray);
+ output->Squeeze();
+
+ for (int xIndex = 0; xIndex < this->ResX; ++xIndex)
+ {
+ delete[] imageBin[xIndex];
+ }
+ delete[] imageBin;
+
+ return 1;
+}
diff --git a/Infovis/Core/vtkContinuousScatterplot.h b/Infovis/Core/vtkContinuousScatterplot.h
new file mode 100644
index 0000000000000000000000000000000000000000..d2c508d0466d105b3eb7ebb3299e29350e370959
--- /dev/null
+++ b/Infovis/Core/vtkContinuousScatterplot.h
@@ -0,0 +1,228 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: vtkPolyDataAlgorithm.h
+
+ 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.
+
+=========================================================================*/
+/**
+* @class vtkContinuousScatterplot
+* @brief Given a 3D domain space represented by an
+* unstructured grid composed of tetrahedral cells with bivariate fields, this filter
+* tessellates each cell in the domain to polyhedral fragments by intersecting the
+* projection of the cell into 2-D range space against two sets of cutting planes, one set
+* is defined along the first field, the second set is defined along the second field. The
+* volume of these subdivided polyhedral fragments can be computed and aggregated over
+* cells to depict the density distribution of the data projection in the bivariate range
+* space.
+*
+* @section Introduction
+* Given a bivariate field (f1,f2) defined on an unstructured grid which is
+* composed of tetrahedral cells, we can initially subdivide each cell based on its
+* projection in the range into a number of fragments along the first field f1, we refer
+* to these polyhedral fragments as Frag(f1) = {frag(f1)_1, frag(f1)_2, ... , frag(f1)_n},
+* where frag(f1)_n refers to the nth fragment along the first field subdivision. Each
+* fragment has a range value and the value difference between the neighbouring fragments
+* is represented as fragment width fw_f1, which is uniformly distributed across the
+* range.
+* Based on the structure of Frag(f1), for each of its cell "frag(f1)_n", we
+* can further subdivide this cell based on the second field f2 using fragment width
+* fw_f2. The tessellation along the second field results in an even finer fragment
+* collection which we refer to as Frag(f1,f2) = {frag(f1,f2)_1, frag(f1,f2)_2, ... ,
+* frag(f1,f2)_m}. We can observe that Frag(f1,f2) is a finer tessellation of the domain
+* than Frag(f1) and will be used to compute the density distribution in the bivariate
+* range space. The algorithm for fragment computation is similar to the first stage of
+* the work in [0].
+* Each fragment "s" in Frag(f1,f2) has range values (f1(s), f2(s)) in the bivariate
+* fields. These values can be further mapped to a 2-D bin with a resolution rexX * resY.
+* The mapped bin index (binIndexX, binIndexY) of the fragment can be computed by linear
+* interpolation on its range values :
+* binIndexX = (int) resX * (f1(s) - f1_min) / (f1_max - f1_min)
+* binIndexY = (int) resY * (f2(s) - f2_min) / (f2_max - f2_min),
+* where (f1_min, f1_max) is the range in first field.
+* Once we know which bin a fragment coincides, the density value in each bin equals to
+* the total geometric volume of the fragments in this bin. This volume distribution
+* over the bins will be exported as a point data array in the output data structure.
+* If we map this 2-D bin to a 2-D image with each bin corresponding to a pixel and
+* bin density to pixel transparency, then the image can be displayed as a continuous
+* scatterplot.
+
+* @section Algorithm
+* The algorithm of this filter can be described as:
+* Require: R.1 The domain space is an unstructured grid data set composed of
+* tetrahedral cells;
+* R.2 The range space contains two scalar fields, say f1 and f2.
+*
+* The most important step is to compute the fragments. The implementation processes
+* the input grid one cell at a time, explicitly computing the intersection of the cell
+* with the cutting planes defined by the fragment boundaries in each scalar field.
+* In order to subdivide the cell, we need to define a list of cutting planes in each
+* field. The interval between neighbouring cutting planes is related to the output 2-D
+* bin resolution (resX, resY) and can be computed as :
+* fw_f1 = (f1_max - f1_min) / resX
+* fw_f2 = (f2_max - f2_min) / resY,
+* where (f1_max,f1_min) is the scalar range of first field.
+*
+* 1. For each tetrahedron T in the input grid:
+*
+* 1.1 Subdivide the cell T based on the first field f1, we will obtain a list
+* of fragments: Frag(f1) = {frag(f1)_1, frag(f1)_2, ... , frag(f1)_n}. The
+* steps for subdivision can be described as:
+*
+* 1.1.1 For each cutting plane s with respect to the first field f1,
+* its field value f1(s) = f1_min + n * fw_f1, where n refers to the n-th
+* cutting plane:
+*
+* 1.1.2. Traverse each edge e starting from point a to b in the cell, we
+* will maintain three data classes, namely fragmentFace,
+* residualFace and cutSet:
+* A. fragmentFace contains vertices in the current fragment.
+* B. cutSet contains vertices whose range values equal to f1(s).
+* This set contains the current cutting plane.
+* C. residualFace contains the rest of the vertices in the cell.
+* In order to classify edge vertices into these classes, the
+* following case table is used for each vertex "a" :
+* case 0 : f1(a)------ f1(s) ------f1(b)
+* condition: f1(a) < f1(s) , f1(b) > f1(s)
+* class: p(s,e), a -> fragmentFace
+* p(s,e) -> cutSet
+* p(s,e) -> residualFace
+*
+* case 1 : f1(b)------ f1(s) ------f1(a)
+* condition: f1(a) > f1(s) , f1(b) < f1(s)
+* class: p(s,e) -> fragmentFace
+* p(s,e) -> cutSet
+* a -> residualFace
+*
+* case 2 : f1(s),f1(a)-------------------f1(b)
+* condition: f1(s) == f1(a), f1(s) <= f1(b)
+* class: a -> fragmentFace
+* a -> residualFace
+* a -> cutSet
+*
+* case 3 : f1(a)-------------------f1(b), f1(s)
+* condition: f1(s) > f1(a), f1(s) == f1(b)
+* class: a -> fragmentFace
+*
+* case 4 : f1(s),f1(b)-------------------f1(a)
+* condition: f1(s) < f1(a), f1(s) == f1(b)
+* class: a -> residualFace
+* Remark: 1. we use "->" to indicate "belongs to" relation.
+* 2. p(s,e) refers to the interpolated point of range value
+* f1(s) on the edge e.
+*
+* 1.1.3. After we have traversed every edge in a cell for the cutting plane
+* s, three classes for storing fragment, cutting plane and residual
+* faces are updated. The faces of the current fragment frag(f1)
+* are the union of all elements in fragmentFace and cutSet.
+*
+* 1.2 Take the output of step 1.1, traverse each fragment in Frag(f1), define a list
+* of cutting planes with respect to field f2, further subdivide the fragments in
+* Frag(f1) following steps from 1.1.2 to 1.1.3. The output of this step will be
+* the fragment collection Frag(f1,f2). Each fragment in Frag(f1,f2) can be further
+* mapped to a 2-D bin based on its range values. The density value in each bin
+* equals to the total geometric volume of the fragments in this bin. This volume
+* distribution over the bins will be exported as a point data array in the output
+* data structure.
+*
+* @section VTK Filter Design
+* The input and output ports of the filter:
+* Input port : the input data set should be a vtkUnstructuredGrid, with each of its
+* cell defined as a tetrahedron. At least two scalar fields are
+* associated with the data. The user needs to specify the name of the
+* two scalar arrays beforehand.
+* Output port: the output data set is a 2D image stored as a vtkImageData.
+* The resolution of the output image can be set by the user.
+* The volume distribution of fragments in each pixel or bin
+* is stored in an point data array named "volume" in the output
+* vtkImageData.
+*
+* @section How To Use This Filter
+* Suppose we have a tetrahedral mesh stored in a vtkUnstructuredGrid, we call this
+* data set "inputData". This data set has two scalar arrays whose names are "f1"
+* and "f2" respectively. We would like the resolution of output image set to (resX,resY).
+* Given these input, this filter can be called as follows in c++ sample code:
+*
+* vtkSmartPointer<vtkContinuousScatterplot> csp =
+* vtkSmartPointer<vtkContinuousScatterplot>::New();
+* csp->SetInputData(inputData);
+* csp->SetField1("f1",resX);
+* csp->SetField2("f2",resY);
+* csp->Update();
+*
+* Then the output, "csp->GetOutput()", will be a vtkImageData containing a scalar
+* array whose name is "volume". This array contains the volume distribution of the
+* fragments.
+*
+* [0] H.Carr and D.Duke, Joint contour nets: Topological analysis of multivariate data.
+* IEEE Transactions on Visualization and Computer Graphics, volume 20,
+* issue 08, pages 1100-1113, 2014
+*/
+
+#ifndef vtkContinuousScatterplot_h
+#define vtkContinuousScatterplot_h
+
+#include "vtkInfovisCoreModule.h" // For export macro
+#include "vtkImageAlgorithm.h"
+
+class VTKINFOVISCORE_EXPORT vtkContinuousScatterplot : public vtkImageAlgorithm
+{
+public:
+ static vtkContinuousScatterplot* New();
+ vtkTypeMacro(vtkContinuousScatterplot, vtkImageAlgorithm);
+ void PrintSelf(ostream& os, vtkIndent indent) VTK_OVERRIDE;
+
+ /**
+ * Get the tolerance used when comparing floating point numbers for equality.
+ */
+ vtkGetMacro(Epsilon, double);
+
+ /**
+ * Set the tolerance used when comparing floating point numbers for equality.
+ */
+ vtkSetMacro(Epsilon, double);
+
+ /**
+ * Specify the name of the first field to be used in subdividing the dataset.
+ * Specify the resolution along x axis of the output image.
+ */
+ void SetField1(char* fieldName, vtkIdType ResX);
+
+ /**
+ * Specify the name of the second field to be used in subdividing the dataset.
+ * Specify the resolution along y axis of the output image.
+ */
+ void SetField2(char* fieldName, vtkIdType ResY);
+
+protected:
+ vtkContinuousScatterplot();
+
+ // Configure input port to accept only vtkUnstructuredGrid.
+ virtual int FillInputPortInformation(int port, vtkInformation* info) VTK_OVERRIDE;
+
+ // Configure out port to be a vtkImageData data set.
+ virtual int FillOutputPortInformation(int port, vtkInformation* info) VTK_OVERRIDE;
+ virtual int RequestData(
+ vtkInformation*, vtkInformationVector**, vtkInformationVector*) VTK_OVERRIDE;
+
+ // Set the tolerance used when comparing floating numbers for equality.
+ double Epsilon;
+
+ // Names of the scalar fields to be used in the filter.
+ char* Fields[2];
+
+ // Resolution of the output image.
+ vtkIdType ResX, ResY;
+
+private:
+ vtkContinuousScatterplot(const vtkContinuousScatterplot&) VTK_DELETE_FUNCTION;
+ void operator=(const vtkContinuousScatterplot&) VTK_DELETE_FUNCTION;
+};
+#endif
diff --git a/Testing/Data/FiberSurface/README.txt.md5 b/Testing/Data/FiberSurface/README.txt.md5
new file mode 100644
index 0000000000000000000000000000000000000000..5374264f5fafeedf20a145530da245769866909b
--- /dev/null
+++ b/Testing/Data/FiberSurface/README.txt.md5
@@ -0,0 +1 @@
+c50c97be6b11d9ee05e7a2f6fde46a0d
diff --git a/Testing/Data/FiberSurface/line_01.vtk.md5 b/Testing/Data/FiberSurface/line_01.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..418951264d24fc18b1eaa7e9f3c8430a6519cecb
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_01.vtk.md5
@@ -0,0 +1 @@
+98fdbe475730fabc6a2934e76b2621cb
diff --git a/Testing/Data/FiberSurface/line_02.vtk.md5 b/Testing/Data/FiberSurface/line_02.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..0c9ca9447d2c38d63d0f66ffb9d31f8fc81e4f09
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_02.vtk.md5
@@ -0,0 +1 @@
+0c890c7ddea677b325c86230629aa69a
diff --git a/Testing/Data/FiberSurface/line_03.vtk.md5 b/Testing/Data/FiberSurface/line_03.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..4a8b6f42900748ea435d8c0b4fc36d29f748e485
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_03.vtk.md5
@@ -0,0 +1 @@
+f7c178c7284c94227a1494fb22d00432
diff --git a/Testing/Data/FiberSurface/line_04.vtk.md5 b/Testing/Data/FiberSurface/line_04.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..67566862844e72dec59dc0fe9b1e970032f6e1ea
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_04.vtk.md5
@@ -0,0 +1 @@
+aae2336f45f5d4742f2e790f0887fbbc
diff --git a/Testing/Data/FiberSurface/line_05.vtk.md5 b/Testing/Data/FiberSurface/line_05.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..17af2690fdfb70d92f1a545dd5676a0761d29df1
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_05.vtk.md5
@@ -0,0 +1 @@
+546fa8373c92e42d28f9fa15e246c393
diff --git a/Testing/Data/FiberSurface/line_11.vtk.md5 b/Testing/Data/FiberSurface/line_11.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..e9f43c2102c5f4588cf95c9c8de2673a8e376727
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_11.vtk.md5
@@ -0,0 +1 @@
+07640d24818c141b3247df9506da9b7d
diff --git a/Testing/Data/FiberSurface/line_12.vtk.md5 b/Testing/Data/FiberSurface/line_12.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..eb03e1f4bf5dba7bdbf352e2cafa6e5382e36024
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_12.vtk.md5
@@ -0,0 +1 @@
+b6edfe4419678ec9d7117bb9f9b44e43
diff --git a/Testing/Data/FiberSurface/line_13.vtk.md5 b/Testing/Data/FiberSurface/line_13.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..b67ff3b551d8f8a7204ae5d064f4def39c2a7e9d
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_13.vtk.md5
@@ -0,0 +1 @@
+e16c7189d8dd0dcac081f3a1b8d72f95
diff --git a/Testing/Data/FiberSurface/line_14.vtk.md5 b/Testing/Data/FiberSurface/line_14.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..86024230aab09c87bf83871960843e624862178e
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_14.vtk.md5
@@ -0,0 +1 @@
+29705575fe5fd9f93441c3f5fed4f4de
diff --git a/Testing/Data/FiberSurface/line_15.vtk.md5 b/Testing/Data/FiberSurface/line_15.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..1d2c330e64baa855de24c7408718e9b74d7e329d
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_15.vtk.md5
@@ -0,0 +1 @@
+445cb1432886b360c2b0ca7d0379f6d0
diff --git a/Testing/Data/FiberSurface/line_21.vtk.md5 b/Testing/Data/FiberSurface/line_21.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..51c7da161dc2bc3b309bbe61ba49c2fa8be55bc0
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_21.vtk.md5
@@ -0,0 +1 @@
+90492c670c1ff8e0fce58b9c810ab235
diff --git a/Testing/Data/FiberSurface/line_22.vtk.md5 b/Testing/Data/FiberSurface/line_22.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..7e086d06a3b0859d02c03ae1362055cc34d40840
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_22.vtk.md5
@@ -0,0 +1 @@
+9c2bbe4f8f0eddf5f7b3ee4e8d711e75
diff --git a/Testing/Data/FiberSurface/line_23.vtk.md5 b/Testing/Data/FiberSurface/line_23.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..c39b71e4889e49b37c9fb3303efc22bdf94db57c
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_23.vtk.md5
@@ -0,0 +1 @@
+4037a6f1f94d2e24a1b80f071c7150b3
diff --git a/Testing/Data/FiberSurface/line_24.vtk.md5 b/Testing/Data/FiberSurface/line_24.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..9034d83b778da4ac00513b4229730eed952b6d8a
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_24.vtk.md5
@@ -0,0 +1 @@
+71527541d8d17d68820ab2ac97c35cdf
diff --git a/Testing/Data/FiberSurface/line_25.vtk.md5 b/Testing/Data/FiberSurface/line_25.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..83bd25725dc718f643008695df064aed86f474a9
--- /dev/null
+++ b/Testing/Data/FiberSurface/line_25.vtk.md5
@@ -0,0 +1 @@
+ab733b35f8938e6e72c1b8697d68d853
diff --git a/Testing/Data/FiberSurface/one_cube.vtk.md5 b/Testing/Data/FiberSurface/one_cube.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..93f2c8f7e7e328dfd1ae7972f74f0d2c147a4455
--- /dev/null
+++ b/Testing/Data/FiberSurface/one_cube.vtk.md5
@@ -0,0 +1 @@
+f5a0f148371af48ecc7d289f0d58342c
diff --git a/Testing/Data/FiberSurface/one_cube_both_forking.vtk.md5 b/Testing/Data/FiberSurface/one_cube_both_forking.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..4660d7b35b6715de4394c0d3e5440c9292fd6852
--- /dev/null
+++ b/Testing/Data/FiberSurface/one_cube_both_forking.vtk.md5
@@ -0,0 +1 @@
+292fbf06bd030d917da38bbae7d1d39d
diff --git a/Testing/Data/FiberSurface/one_cube_closed.vtk.md5 b/Testing/Data/FiberSurface/one_cube_closed.vtk.md5
new file mode 100644
index 0000000000000000000000000000000000000000..f81705f58902913dab406e3a03520635c45405f7
--- /dev/null
+++ b/Testing/Data/FiberSurface/one_cube_closed.vtk.md5
@@ -0,0 +1 @@
+460894a81078cef1eb159e1dbb24bc01
diff --git a/Testing/Data/cube.vtu.md5 b/Testing/Data/cube.vtu.md5
new file mode 100644
index 0000000000000000000000000000000000000000..b73b319e6e4bafd355d16e5f972671601bba0ce2
--- /dev/null
+++ b/Testing/Data/cube.vtu.md5
@@ -0,0 +1 @@
+876bdbeba77529faf39aca1503b8d201