Commit from Revision 1 of VTK Journal article

http://www.midasjournal.org/browse/publication/792

CMakeLists.txt

@@ -0,0 +1,59 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
+INCLUDE(${CMAKE_SOURCE_DIR}/IJMacros.txt)
+
+#Change PROJECT_NAME to the name of your project
+PROJECT(PartialVolumeModeller)
+
+#The following lines are required to use Dart
+ENABLE_TESTING()
+INCLUDE(Dart)
+
+#Declare any external dependencies that your project may have here.
+#examples include: ITK, VTK, JPEG, PNG, OpenGL, ZLIB, Perl, Java
+#If you're not sure what name to use, look in the Modules directory of your
+#cmake install and check that a file named Find(Package).cmake exists
+#
+# The packages can be specified with a version number, for example:
+#
+# ITK 2.8.1
+# ITK 3.2.0
+#
+# If no version is specified, the most recent release of the package
+# will be used.
+SET(Required_Packages
+  VTK
+)
+
+#this foreach loads all of the packages that you specified as required.
+#It shouldn't need to be modified.
+FOREACH(Package ${Required_Packages})
+  LOADPACKAGE(${Package})
+ENDFOREACH(Package)
+
+#Set any libraries that your project depends on.
+#examples: ITKCommon, VTKRendering, etc
+SET(Libraries
+  vtkFiltering
+  vtkGraphics
+  vtkImaging
+  vtkIO
+)
+
+#the following block of code is an example of how to build an executable in
+#cmake.  Unmodified, it will add an executable called "MyExe" to the project.
+#MyExe will be built using the files MyClass.h and MyClass.cxx, and it will
+#be linked to all the libraries you specified above.
+#You can build more than one executable per project
+SET(CurrentExe "TestPartialVolumeModeller")
+ADD_EXECUTABLE(${CurrentExe}
+  vtkPartialVolumeModeller.cxx
+  TestPartialVolumeModeller.cxx)
+TARGET_LINK_LIBRARIES(${CurrentExe} ${Libraries})
+
+#the following line is an example of how to add a test to your project.
+#Testname is the title for this particular test.  ExecutableToRun is the
+#program which will be running this test.  It can either be a part of this
+#project or an external executable.  After that list any args that are needed
+#for this test.  Include as many tests as you like.  If your project doesn't have
+#any tests you can comment out or delete the following line.
+ADD_TEST(TestPartialVolumeModeller ${CurrentExe})

IJMacros.txt

@@ -0,0 +1,79 @@
+#Macro to find a package and load it
+#(you shouldn't need to modify this code)
+MACRO(LOADPACKAGE Package)
+  SET(Included FALSE)
+  IF(EXISTS "/home/tester/XMLTestParser.py")
+    #If we're in the test environment, check and see if we're being asked for
+    #a specific version of some package.  If so, we'll provide a direct path
+    #instead of counting on CMake to choose the right version.
+    IF(${Package} MATCHES "^ITK.*3[.]2[.]0$")
+      SET(ITK_DIR "/home/tester/ITK3.2.0/bin")
+      INCLUDE("/usr/local/share/CMake/Modules/FindITK.cmake")
+      INCLUDE(${ITK_USE_FILE})
+      SET(Included TRUE)
+    ENDIF(${Package} MATCHES "^ITK.*3[.]2[.]0$")
+    IF(${Package} MATCHES "^ITK.*2[.]2[.]1$")
+      SET(ITK_DIR "/home/tester/ITK2.2.1/bin")
+      INCLUDE("/usr/local/share/CMake/Modules/FindITK.cmake")
+      INCLUDE(${ITK_USE_FILE})
+      SET(Included TRUE)
+    ENDIF(${Package} MATCHES "^ITK.*2[.]2[.]1$")
+    IF(NOT Included AND ${Package} MATCHES "^ITK.*1[.]8[.]1$")
+      SET(ITK_DIR "/home/tester/ITK1.8.1/bin")
+      INCLUDE("/usr/local/share/CMake/Modules/FindITK.cmake")
+      INCLUDE(${ITK_USE_FILE})
+      SET(Included TRUE)
+    ENDIF(NOT Included AND ${Package} MATCHES "^ITK.*1[.]8[.]1$")
+    IF(NOT Included AND ${Package} MATCHES "^VTK.*4[.]4$")
+      SET(VTK_DIR "/home/tester/VTK4.4/bin")
+      INCLUDE("/usr/local/share/CMake/Modules/FindVTK.cmake")
+      INCLUDE(${VTK_USE_FILE})
+      SET(Included TRUE)
+    ENDIF(NOT Included AND ${Package} MATCHES "^VTK.*4[.]4$")
+    IF(NOT Included AND ${Package} MATCHES "^VTK.*5[.]0$")
+      SET(VTK_DIR "/home/tester/VTK5.0/bin")
+      INCLUDE("/usr/local/share/CMake/Modules/FindVTK.cmake")
+      INCLUDE(${VTK_USE_FILE})
+      SET(Included TRUE)
+    ENDIF(NOT Included AND ${Package} MATCHES "^VTK.*5[.]0$")
+    #If we get this far and we still haven't found a match, set it up so the
+    #next block of code has a chance at finding the package.
+    IF(NOT Included AND ${Package} MATCHES "^VTK")
+      SET(Package "VTK")
+    ENDIF(NOT Included AND ${Package} MATCHES "^VTK")
+    IF(NOT Included AND ${Package} MATCHES "^ITK")
+      SET(Package "ITK")
+    ENDIF(NOT Included AND ${Package} MATCHES "^ITK")
+  ENDIF(EXISTS "/home/tester/XMLTestParser.py")
+
+  #no point in executing the code below if we already found the package we're
+  #looking for.
+  IF(NOT Included)
+    FIND_PACKAGE(${Package})
+    IF(${Package}_FOUND)
+      #most packages define a Package_INCLUDE_DIR variable, so we'll check for
+      #that first
+      IF(${Package}_INCLUDE_DIR)
+  INCLUDE(${${Package}_INCLUDE_DIR})
+  SET(Included TRUE)
+      ELSE(${Package}_INCLUDE_DIR)
+  #VTK and ITK prefer to define a Package_USE_FILE, so we need to look for
+  #that too
+  IF(${Package}_USE_FILE)
+    INCLUDE(${${Package}_USE_FILE})
+    SET(Included TRUE)
+  ENDIF(${Package}_USE_FILE)
+      ENDIF(${Package}_INCLUDE_DIR)
+      #then there's some other pesky packages that don't like to define standard
+      #variables at all.  If you're trying to include one of those you might have
+      #to do a little bit of investigating on your own.
+      IF(NOT Included)
+  MESSAGE(FATAL_ERROR "${Package} was found, but couldn't be included.\n
+    Try including it manually out of the FOREACH in the CMakeLists file.\n
+    Look at Find${Package}.cmake in the CMake module directory for clues
+    on what you're supposed to include.  Good luck.")
+      ENDIF(NOT Included)
+    ENDIF(${Package}_FOUND)
+  ENDIF(NOT Included)
+ENDMACRO(LOADPACKAGE)
+

TestPartialVolumeModeller.cxx

@@ -0,0 +1,179 @@
+/*=========================================================================
+
+  Program:   Visualization Toolkit
+  Module:    vtkPartialVolumeModeller.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 <iostream>
+
+#include <vtkDoubleArray.h>
+#include <vtkImageData.h>
+#include <vtkImageShiftScale.h>
+#include <vtkPNGWriter.h>
+#include <vtkPoints.h>
+#include <vtkRectilinearGridToTetrahedra.h>
+#include <vtkTetra.h>
+#include <vtkType.h>
+#include <vtkUnstructuredGrid.h>
+#include <vtkVoxelModeller.h>
+#include <vtkXMLImageDataWriter.h>
+#include <vtkXMLUnstructuredGridWriter.h>
+
+#include "vtkPartialVolumeModeller.h"
+
+
+int main(int argc, char* argv[])
+{
+  // Create a box partitioned into tetrahedra.
+  double spacing = 5.0;
+  double boxVoxelsCoveredX = 20;
+  double boxVoxelsCoveredY = 30;
+  double boxVoxelsCoveredZ = 15;
+
+  double boxWidth = spacing * boxVoxelsCoveredX;
+  double boxHeight = spacing * boxVoxelsCoveredY;
+  double boxDepth = spacing * boxVoxelsCoveredZ;
+
+  // Node numbering scheme
+  //       7-----6
+  //      /|    /|
+  //     4-+---5 |
+  //     | 3---+-2
+  //     |/    |/
+  //     0-----1
+  double tetPoints[8][3] = {
+    {0.0, 0.0, boxDepth},
+    {boxWidth, 0.0, boxDepth},
+    {boxWidth, 0.0, 0.0},
+    {0.0, 0.0, 0.0},
+    {0.0, boxHeight, boxDepth},
+    {boxWidth, boxHeight, boxDepth},
+    {boxWidth, boxHeight, 0.0},
+    {0.0, boxHeight, 0.0}};
+
+  vtkIdType tetIndices[6][4] = {
+    {1, 3, 0, 4},
+    {4, 1, 3, 5},
+    {4, 7, 5, 3},
+    {5, 7, 6, 3},
+    {1, 3, 5, 2},
+    {3, 2, 6, 5}
+  };
+
+  vtkUnstructuredGrid *grid = vtkUnstructuredGrid::New();
+
+  vtkPoints* gridPoints = vtkPoints::New();
+
+  for (int i = 0; i < 8; i++)
+    {
+    gridPoints->InsertNextPoint(tetPoints[i]);
+    }
+  grid->SetPoints(gridPoints);
+
+  for (int i = 0; i < 6; i++)
+    {
+    grid->InsertNextCell(VTK_TETRA, 4, tetIndices[i]);
+    }
+
+  vtkXMLUnstructuredGridWriter *gridWriter = vtkXMLUnstructuredGridWriter::New();
+  gridWriter->SetFileName("TetrahedralGrid.vtu");
+  gridWriter->SetInput(grid);
+  gridWriter->Write();
+
+  // Set up partial volume modeller to create an image larger than the
+  // extent of the box.
+  double padding = 3 * spacing;
+  double off = 0.0; // An offset applied to all modeller bounds
+  double xmin = -padding + off;
+  double xmax = boxWidth + padding + off;
+  double ymin = -padding + off;
+  double ymax = boxHeight + padding + off;
+  double zmin = -padding + off;
+  double zmax = boxDepth + padding + off;
+
+  int imageDimsX = static_cast<int>((xmax - xmin) / spacing) + 1;
+  int imageDimsY = static_cast<int>((ymax - ymin) / spacing) + 1;
+  int imageDimsZ = static_cast<int>((zmax - zmin) / spacing) + 1;
+
+  // Try out the vtkPartialVolumeModeller
+  vtkPartialVolumeModeller *partialVolumeModeller = vtkPartialVolumeModeller::New();
+  partialVolumeModeller->SetModelBounds(xmin, xmax, ymin, ymax, zmin, zmax);
+  partialVolumeModeller->SetSampleDimensions(imageDimsX, imageDimsY, imageDimsZ);
+  partialVolumeModeller->SetInput(grid);
+  partialVolumeModeller->Update();
+
+  // Check out the results at the side edge of the box. Should be 0.25.
+  double edgeValue = partialVolumeModeller->GetOutput()->
+    GetScalarComponentAsDouble(3, 3, 4, 0);
+  if (fabs(edgeValue - 0.25) > 1e-5)
+    {
+    std::cerr << "Expected value at voxel (3, 3, 4) should be 0.25, got "
+              << edgeValue << std::endl;
+    return EXIT_FAILURE;
+    }
+
+  // Check out the results at the corner of the box. Should be 0.125.
+  double cornerValue = partialVolumeModeller->GetOutput()->
+    GetScalarComponentAsDouble(3, 3, 3, 0);
+  if (fabs(cornerValue - 0.125) > 1e-5)
+    {
+    std::cerr << "Expected value at voxel (3, 3, 3) to be 0.125, got "
+              << cornerValue << std::endl;
+    return EXIT_FAILURE;
+    }
+
+#if 0
+  vtkXMLImageDataWriter* imageWriter = vtkXMLImageDataWriter::New();
+  imageWriter->SetFileName("PartialVolumeImageData.vti");
+  imageWriter->SetInputConnection(partialVolumeModeller->GetOutputPort());
+  imageWriter->Write();
+
+  // Scale values to range necessary for image file formats.
+  vtkImageShiftScale *shiftScale = vtkImageShiftScale::New();
+  shiftScale->SetShift(0.0);
+  shiftScale->SetScale(255.0);
+  shiftScale->SetOutputScalarTypeToUnsignedChar();
+  shiftScale->SetInputConnection(partialVolumeModeller->GetOutputPort());
+
+  // Save as a PNG.
+  vtkPNGWriter *writer = vtkPNGWriter::New();
+  writer->SetFilePattern("PartialVolumeModellerImage%04d.png");
+  writer->SetFileDimensionality(2);
+  writer->SetInputConnection(shiftScale->GetOutputPort());
+  writer->Write();
+
+  // Try the vtkVoxelModeller for comparison
+  vtkVoxelModeller *voxelModeller = vtkVoxelModeller::New();
+  voxelModeller->SetModelBounds(xmin, xmax, ymin, ymax, zmin, zmax);
+  voxelModeller->SetSampleDimensions(imageDimsX, imageDimsY, imageDimsZ);
+  voxelModeller->SetScalarTypeToFloat();
+  voxelModeller->SetInput(grid);
+
+  imageWriter->SetFileName("VoxelModellerImageData.vti");
+  imageWriter->SetInputConnection(partialVolumeModeller->GetOutputPort());
+  imageWriter->Write();
+
+  shiftScale->SetInputConnection(voxelModeller->GetOutputPort());
+
+  writer->SetFilePattern("VoxelModellerImage%04d.png");
+  writer->Write();
+
+  writer->Delete();
+  imageWriter->Delete();
+  shiftScale->Delete();
+  voxelModeller->Delete();
+#endif
+
+  partialVolumeModeller->Delete();
+
+  return EXIT_SUCCESS;
+}