From 30100e2ae83e349f2843de50b2c2354a7f722225 Mon Sep 17 00:00:00 2001 From: Christopher Meyer Sewell - 188584 Date: Tue, 1 Sep 2015 17:09:36 -0600 Subject: [PATCH] Adding examples directory with isosurface rendering example --- CMakeLists.txt | 7 + examples/CMakeLists.txt | 44 ++++ examples/IsosurfaceUniformGrid.cxx | 249 ++++++++++++++++++ examples/quaternion.h | 79 ++++++ vtkm/VectorAnalysis.h | 1 + vtkm/worklet/IsosurfaceUniformGrid.h | 18 +- .../testing/UnitTestIsosurfaceUniformGrid.cxx | 3 +- 7 files changed, 398 insertions(+), 3 deletions(-) create mode 100644 examples/CMakeLists.txt create mode 100644 examples/IsosurfaceUniformGrid.cxx create mode 100755 examples/quaternion.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 4da2f8e10..179dbda8c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -70,6 +70,7 @@ option(VTKm_ENABLE_TESTING "Enable VTKm Testing" ON) option(VTKm_ENABLE_BENCHMARKS "Enable VTKm Benchmarking" OFF) option(VTKm_BUILD_DOCUMENTATION "Build Doxygen documentation" OFF) +option(VTKm_BUILD_EXAMPLES "Build examples" OFF) option(VTKm_USE_DOUBLE_PRECISION "Use double precision for floating point calculations" @@ -194,6 +195,12 @@ if (VTKm_BUILD_DOCUMENTATION) vtkm_build_documentation() endif() +#----------------------------------------------------------------------------- +# Build examples +if(VTKm_BUILD_EXAMPLES) + add_subdirectory(examples) +endif(VTKm_BUILD_EXAMPLES) + #----------------------------------------------------------------------------- # Configuration for build directory. set(VTKm_INCLUDE_DIRS_CONFIG "${VTKm_SOURCE_DIR};${VTKm_BINARY_DIR}") diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 000000000..1b5249665 --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,44 @@ +##============================================================================= +## +## Copyright (c) Kitware, Inc. +## All rights reserved. +## See LICENSE.txt for details. +## +## This software is distributed WITHOUT ANY WARRANTY; without even +## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +## PURPOSE. See the above copyright notice for more information. +## +## Copyright 2015 Sandia Corporation. +## Copyright 2015 UT-Battelle, LLC. +## Copyright 2015 Los Alamos National Security. +## +## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +## the U.S. Government retains certain rights in this software. +## Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +## Laboratory (LANL), the U.S. Government retains certain rights in +## this software. +## +##============================================================================= + +cmake_minimum_required(VERSION 2.8.11) +find_package(OpenGL) +find_package(GLUT) +find_package(TBB) + +if(OPENGL_FOUND AND GLUT_FOUND) + add_executable(IsosurfaceUniformGrid_SERIAL IsosurfaceUniformGrid.cxx) + target_include_directories(IsosurfaceUniformGrid_SERIAL PRIVATE ${GLUT_INCLUDE_DIR} ${OPENGL_INCLUDE_DIR}) + target_link_libraries(IsosurfaceUniformGrid_SERIAL ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES}) + + if(VTKm_Cuda_FOUND) + cuda_add_executable(IsosurfaceUniformGrid_CUDA IsosurfaceUniformGrid.cxx) + target_link_libraries(IsosurfaceUniformGrid_CUDA ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES}) + endif() + + if(VTKm_ENABLE_TBB) + add_executable(IsosurfaceUniformGrid_TBB IsosurfaceUniformGrid.cxx) + target_include_directories(IsosurfaceUniformGrid_TBB PRIVATE ${GLUT_INCLUDE_DIR} ${OPENGL_INCLUDE_DIR} ${TBB_INCLUDE_DIRS}) + target_link_libraries(IsosurfaceUniformGrid_TBB ${OPENGL_LIBRARIES} ${GLUT_LIBRARIES} ${TBB_LIBRARIES}) + endif() + +endif() diff --git a/examples/IsosurfaceUniformGrid.cxx b/examples/IsosurfaceUniformGrid.cxx new file mode 100644 index 000000000..9f125e3e2 --- /dev/null +++ b/examples/IsosurfaceUniformGrid.cxx @@ -0,0 +1,249 @@ +//============================================================================ +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2014 Sandia Corporation. +// Copyright 2014 UT-Battelle, LLC. +// Copyright 2014 Los Alamos National Security. +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +//============================================================================ + +#include +#include + +#include +#include + +#include + +#if defined (__APPLE__) +# include +#else +# include +#endif + +#include "quaternion.h" + +#include + +typedef VTKM_DEFAULT_DEVICE_ADAPTER_TAG DeviceAdapter; + +vtkm::worklet::IsosurfaceFilterUniformGrid *isosurfaceFilter; +vtkm::cont::ArrayHandle > verticesArray, normalsArray; +vtkm::cont::ArrayHandle scalarsArray; +Quaternion qrot; +int lastx, lasty; +int mouse_state = 1; + +namespace { + +class TangleField : public vtkm::worklet::WorkletMapField +{ +public: + typedef void ControlSignature(FieldIn vertexId, FieldOut v); + typedef void ExecutionSignature(_1, _2); + typedef _1 InputDomain; + + const vtkm::Id xdim, ydim, zdim; + const float xmin, ymin, zmin, xmax, ymax, zmax; + const vtkm::Id cellsPerLayer; + + VTKM_CONT_EXPORT + TangleField(const vtkm::Id3 dims, const float mins[3], const float maxs[3]) : xdim(dims[0]), ydim(dims[1]), zdim(dims[2]), + xmin(mins[0]), ymin(mins[1]), zmin(mins[2]), xmax(maxs[0]), ymax(maxs[1]), zmax(maxs[2]), cellsPerLayer((xdim) * (ydim)) { }; + + VTKM_EXEC_EXPORT + void operator()(const vtkm::Id &vertexId, vtkm::Float32 &v) const + { + const vtkm::Id x = vertexId % (xdim); + const vtkm::Id y = (vertexId / (xdim)) % (ydim); + const vtkm::Id z = vertexId / cellsPerLayer; + + const float fx = static_cast(x) / static_cast(xdim-1); + const float fy = static_cast(y) / static_cast(xdim-1); + const float fz = static_cast(z) / static_cast(xdim-1); + + const vtkm::Float32 xx = 3.0f*(xmin+(xmax-xmin)*(fx)); + const vtkm::Float32 yy = 3.0f*(ymin+(ymax-ymin)*(fy)); + const vtkm::Float32 zz = 3.0f*(zmin+(zmax-zmin)*(fz)); + + v = (xx*xx*xx*xx - 5.0f*xx*xx + yy*yy*yy*yy - 5.0f*yy*yy + zz*zz*zz*zz - 5.0f*zz*zz + 11.8f) * 0.2f + 0.5f; + } +}; + + +vtkm::cont::DataSet MakeIsosurfaceTestDataSet(vtkm::Id3 dims) +{ + vtkm::cont::DataSet dataSet; + + const vtkm::Id3 vdims(dims[0] + 1, dims[1] + 1, dims[2] + 1); + const vtkm::Id dim3 = dims[0]*dims[1]*dims[2]; + + float mins[3] = {-1.0f, -1.0f, -1.0f}; + float maxs[3] = {1.0f, 1.0f, 1.0f}; + + vtkm::cont::ArrayHandle fieldArray; + vtkm::cont::ArrayHandleCounting vertexCountImplicitArray(0, vdims[0]*vdims[1]*vdims[2]); + vtkm::worklet::DispatcherMapField tangleFieldDispatcher(TangleField(vdims, mins, maxs)); + tangleFieldDispatcher.Invoke(vertexCountImplicitArray, fieldArray); + + vtkm::cont::ArrayHandleUniformPointCoordinates coordinates(vdims); + dataSet.AddCoordinateSystem( + vtkm::cont::CoordinateSystem("coordinates", 1, coordinates)); + + dataSet.AddField(vtkm::cont::Field("nodevar", 1, vtkm::cont::Field::ASSOC_POINTS, fieldArray)); + + std::vector cellvar( static_cast(dim3) ); + for(std::size_t i=0; i < cellvar.size(); i++) + { + cellvar[i] = vtkm::Float32(i); + } + + vtkm::cont::Field cellField("cellvar", 1, + vtkm::cont::Field::ASSOC_CELL_SET, + "cells", + cellvar); + dataSet.AddField(cellField); + + static const vtkm::IdComponent ndim = 3; + vtkm::cont::CellSetStructured cellSet("cells"); + cellSet.SetPointDimensions(vdims); + dataSet.AddCellSet(cellSet); + + return dataSet; +} + +} + + +void initializeGL() +{ + glClearColor(0.0f, 0.0f, 0.0f, 0.0f); + glEnable(GL_DEPTH_TEST); + glShadeModel(GL_SMOOTH); + + float white[] = { 0.8, 0.8, 0.8, 1.0 }; + float black[] = { 0.0, 0.0, 0.0, 1.0 }; + float lightPos[] = { 10.0, 10.0, 10.5, 1.0 }; + + glLightfv(GL_LIGHT0, GL_AMBIENT, white); + glLightfv(GL_LIGHT0, GL_DIFFUSE, white); + glLightfv(GL_LIGHT0, GL_SPECULAR, black); + glLightfv(GL_LIGHT0, GL_POSITION, lightPos); + + glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 1); + + glEnable(GL_LIGHTING); + glEnable(GL_LIGHT0); + glEnable(GL_NORMALIZE); + glColorMaterial(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE); + glEnable(GL_COLOR_MATERIAL); +} + + +void displayCall() +{ + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glEnable(GL_DEPTH_TEST); + + glMatrixMode(GL_PROJECTION); + glLoadIdentity(); + gluPerspective( 45.0f, 1.0f, 1.0f, 20.0f); + + glMatrixMode(GL_MODELVIEW); + glLoadIdentity(); + gluLookAt(0.0f, 0.0f, 3.0f, 0.0f, 0.0f, 0.0f, 0.0f, 1.0f, 0.0f); + + glPushMatrix(); + float rotationMatrix[16]; + qrot.getRotMat(rotationMatrix); + glMultMatrixf(rotationMatrix); + + glColor3f(0.1f, 0.1f, 0.6f); + + glBegin(GL_TRIANGLES); + for (unsigned int i=0; i curNormal = normalsArray.GetPortalConstControl().Get(i); + vtkm::Vec curVertex = verticesArray.GetPortalConstControl().Get(i); + glNormal3f(curNormal[0], curNormal[1], curNormal[2]); + glVertex3f(curVertex[0], curVertex[1], curVertex[2]); + } + glEnd(); + + glPopMatrix(); + glutSwapBuffers(); +} + + +void mouseMove(int x, int y) +{ + int dx = x - lastx; + int dy = y - lasty; + + if (mouse_state == 0) + { + Quaternion newRotX; + newRotX.setEulerAngles(-0.2*dx*M_PI/180.0, 0.0, 0.0); + qrot.mul(newRotX); + + Quaternion newRotY; + newRotY.setEulerAngles(0.0, 0.0, -0.2*dy*M_PI/180.0); + qrot.mul(newRotY); + } + lastx = x; + lasty = y; + + glutPostRedisplay(); +} + + +void mouseCall(int button, int state, int x, int y) +{ + if (button == 0) mouse_state = state; + if ((button == 0) && (state == 0)) { lastx = x; lasty = y; } +} + + +int main(int argc, char* argv[]) +{ + std::cout << "IsosurfaceUniformGrid Example" << std::endl; + + vtkm::Id3 dims(16,16,16); + vtkm::cont::DataSet dataSet = MakeIsosurfaceTestDataSet(dims); + + isosurfaceFilter = new vtkm::worklet::IsosurfaceFilterUniformGrid(dims, dataSet); + + isosurfaceFilter->Run(0.5, + dataSet.GetField("nodevar").GetData(), + verticesArray, + normalsArray, + scalarsArray); + + std::cout << verticesArray.GetNumberOfValues() << std::endl; + + lastx = lasty = 0; + glutInit(&argc, argv); + glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH); + glutInitWindowSize(1000, 1000); + glutInitWindowPosition(300, 200); + glutCreateWindow("VTK-m Isosurface"); + initializeGL(); + glutDisplayFunc(displayCall); + glutMotionFunc(mouseMove); + glutMouseFunc(mouseCall); + glutMainLoop(); + + return 0; +} + diff --git a/examples/quaternion.h b/examples/quaternion.h new file mode 100755 index 000000000..c2cdd7105 --- /dev/null +++ b/examples/quaternion.h @@ -0,0 +1,79 @@ +//============================================================================= +// +// Copyright (c) Kitware, Inc. +// All rights reserved. +// See LICENSE.txt for details. +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notice for more information. +// +// Copyright 2015 Sandia Corporation. +// Copyright 2015 UT-Battelle, LLC. +// Copyright 2015 Los Alamos National Security. +// +// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, +// the U.S. Government retains certain rights in this software. +// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National +// Laboratory (LANL), the U.S. Government retains certain rights in +// this software. +// +//============================================================================= + +/* + * quaternion.h + * + * Created on: Oct 10, 2014 + * Author: csewell + */ + +#ifndef QUATERNION_H_ +#define QUATERNION_H_ + +#include +#include + +class Quaternion +{ +public: + Quaternion() { x = y = z = 0.0; w = 1.0; } + Quaternion(double x, double y, double z, double w) : x(x), y(y), z(z), w(w) {}; + void set(double ax, double ay, double az, double aw) { x = ax; y = ay; z = az; w = aw; } + void normalize() + { + float norm = sqrt(x*x + y*y + z*z + w*w); + if (norm > 0.00001) { x /= norm; y /= norm; z /= norm; w /= norm; } + } + void mul(Quaternion q) + { + float tx, ty, tz, tw; + tx = w*q.x + x*q.w + y*q.z - z*q.y; + ty = w*q.y + y*q.w + z*q.x - x*q.z; + tz = w*q.z + z*q.w + x*q.y - y*q.x; + tw = w*q.w - x*q.x - y*q.y - z*q.z; + + x = tx; y = ty; z = tz; w = tw; + } + void setEulerAngles(float pitch, float yaw, float roll) + { + w = cos(pitch/2.0)*cos(yaw/2.0)*cos(roll/2.0) - sin(pitch/2.0)*sin(yaw/2.0)*sin(roll/2.0); + x = sin(pitch/2.0)*sin(yaw/2.0)*cos(roll/2.0) + cos(pitch/2.0)*cos(yaw/2.0)*sin(roll/2.0); + y = sin(pitch/2.0)*cos(yaw/2.0)*cos(roll/2.0) + cos(pitch/2.0)*sin(yaw/2.0)*sin(roll/2.0); + z = cos(pitch/2.0)*sin(yaw/2.0)*cos(roll/2.0) - sin(pitch/2.0)*cos(yaw/2.0)*sin(roll/2.0); + + normalize(); + } + + void getRotMat(float* m) const + { + for (int i=0; i<16; i++) m[i] = 0.0; m[15] = 1.0; + m[0] = 1 - 2*y*y - 2*z*z; m[1] = 2*x*y - 2*z*w; m[2] = 2*x*z + 2*y*w; + m[4] = 2*x*y + 2*z*w; m[5] = 1 - 2*x*x - 2*z*z; m[6] = 2*y*z - 2*x*w; + m[8] = 2*x*z - 2*y*w; m[9] = 2*y*z + 2*x*w; m[10] = 1 - 2*x*x - 2*y*y; + } + + double x,y,z,w; +}; + + +#endif /* QUATERNION_H_ */ diff --git a/vtkm/VectorAnalysis.h b/vtkm/VectorAnalysis.h index 03738931e..6c21f01de 100644 --- a/vtkm/VectorAnalysis.h +++ b/vtkm/VectorAnalysis.h @@ -139,6 +139,7 @@ RMagnitudeTemplate(const T &x, vtkm::TypeTraitsVectorTag) /// as fast as MagnitudeSquared. /// template +VTKM_EXEC_CONT_EXPORT typename vtkm::VecTraits::ComponentType RMagnitude(const T &x) { diff --git a/vtkm/worklet/IsosurfaceUniformGrid.h b/vtkm/worklet/IsosurfaceUniformGrid.h index 6ade4ed94..ca7ea7cd0 100644 --- a/vtkm/worklet/IsosurfaceUniformGrid.h +++ b/vtkm/worklet/IsosurfaceUniformGrid.h @@ -105,6 +105,7 @@ public: typedef typename vtkm::cont::ArrayHandle >::template ExecutionTypes::Portal VectorPortalType; VectorPortalType Vertices; + VectorPortalType Normals; typedef typename vtkm::cont::ArrayHandle IdArrayHandle; typedef typename IdArrayHandle::ExecutionTypes::PortalConst IdPortalType; @@ -115,14 +116,15 @@ public: template VTKM_CONT_EXPORT IsoSurfaceGenerate(const float ivalue, const vtkm::Id3 cdims, IdPortalType triTablePortal, - const U & field, const U & source, const W & vertices, const X & scalars) : + const U & field, const U & source, const W & vertices, const W & normals, const X & scalars) : Isovalue(ivalue), xdim(cdims[0]), ydim(cdims[1]), zdim(cdims[2]), xmin(-1), ymin(-1), zmin(-1), xmax(1), ymax(1), zmax(1), Field( field.PrepareForInput( DeviceAdapter() ) ), Source( source.PrepareForInput( DeviceAdapter() ) ), - Scalars(scalars), Vertices(vertices), + Normals(normals), + Scalars(scalars), TriTable(triTablePortal), cellsPerLayer(xdim * ydim), pointsPerLayer ((xdim+1)*(ydim+1)) @@ -224,6 +226,16 @@ public: this->Vertices.Set(outputVertId + v, vtkm::Lerp(p[v0], p[v1], t)); this->Scalars.Set(outputVertId + v, vtkm::Lerp(s[v0], s[v1], t)); } + + vtkm::Vec vertex0 = this->Vertices.Get(outputVertId + 0); + vtkm::Vec vertex1 = this->Vertices.Get(outputVertId + 1); + vtkm::Vec vertex2 = this->Vertices.Get(outputVertId + 2); + + vtkm::Vec curNorm = vtkm::Cross(vertex1-vertex0, vertex2-vertex0); + vtkm::Normalize(curNorm); + this->Normals.Set(outputVertId + 0, curNorm); + this->Normals.Set(outputVertId + 1, curNorm); + this->Normals.Set(outputVertId + 2, curNorm); } }; @@ -242,6 +254,7 @@ public: void Run(const float &isovalue, IsoField isoField, vtkm::cont::ArrayHandle< vtkm::Vec > &verticesArray, + vtkm::cont::ArrayHandle< vtkm::Vec > &normalsArray, vtkm::cont::ArrayHandle &scalarsArray) { typedef typename vtkm::cont::DeviceAdapterAlgorithm DeviceAlgorithms; @@ -303,6 +316,7 @@ public: field, field, verticesArray.PrepareForOutput(numTotalVertices, DeviceAdapter()), + normalsArray.PrepareForOutput(numTotalVertices, DeviceAdapter()), scalarsArray.PrepareForOutput(numTotalVertices, DeviceAdapter()) ); diff --git a/vtkm/worklet/testing/UnitTestIsosurfaceUniformGrid.cxx b/vtkm/worklet/testing/UnitTestIsosurfaceUniformGrid.cxx index 1d5944f84..1aaea9590 100644 --- a/vtkm/worklet/testing/UnitTestIsosurfaceUniformGrid.cxx +++ b/vtkm/worklet/testing/UnitTestIsosurfaceUniformGrid.cxx @@ -121,11 +121,12 @@ void TestIsosurfaceUniformGrid() vtkm::worklet::IsosurfaceFilterUniformGrid isosurfaceFilter(dims, dataSet); - vtkm::cont::ArrayHandle > verticesArray; + vtkm::cont::ArrayHandle > verticesArray, normalsArray; vtkm::cont::ArrayHandle scalarsArray; isosurfaceFilter.Run(0.5, dataSet.GetField("nodevar").GetData(), verticesArray, + normalsArray, scalarsArray); VTKM_TEST_ASSERT(test_equal(verticesArray.GetNumberOfValues(), 480),