Apply master bug fixes and dependency updates to 6.8

.github/workflows/download_test_data.sh

@@ -18,17 +18,19 @@ if [ "$PLATFORM" == "windows" ]; then
     choco install wget --no-progress
 fi
 
-# Tiny test data downloads
+# Tiny test data. A failed attempt leaves a partial file behind, so retry on whether the
+# archive is complete rather than on whether it exists, and let wget resume it.
+cd /tmp
 for run in {1..5}; do
-    if [ ! -f /tmp/tiny_test_cache.tar.gz ] ; then
-	cd /tmp
-	wget https://www.sci.utah.edu/~shapeworks/doc-resources/zips/tiny_test_cache.tar.gz || true
-	sleep 1
+    if tar -tzf /tmp/tiny_test_cache.tar.gz > /dev/null 2>&1 ; then
+	break
     fi
+    wget --continue --tries=3 --timeout=60 https://www.sci.utah.edu/~shapeworks/doc-resources/zips/tiny_test_cache.tar.gz || true
+    sleep 1
 done
 
-# Decompress if successful
-if [ -f /tmp/tiny_test_cache.tar.gz ] ; then
+# Only extract a complete archive: a half-extracted cache sends the use cases to the web.
+if tar -tzf /tmp/tiny_test_cache.tar.gz > /dev/null 2>&1 ; then
     cd $BASE/Examples/Python
     tar -xzvf /tmp/tiny_test_cache.tar.gz
     rm /tmp/tiny_test_cache.tar.gz

Applications/shapeworks/shapeworks.cpp

@@ -6,11 +6,14 @@
 #include "Executable.h"
 #include "Commands.h"
 #include <Profiling.h>
+#include <ShapeworksUtils.h>
 
 using namespace shapeworks;
 
 int main(int argc, char *argv[])
 {
+  ShapeWorksUtils::initialize_path_handling();
+
   Executable shapeworks;
 
   // Image Commands

Libs/Analyze/Shape.cpp

@@ -5,6 +5,7 @@
 #include <Particles/ParticleSystemEvaluation.h>
 #include <Project/ProjectUtils.h>
 #include <Shape.h>
+#include <ShapeworksUtils.h>
 #include <Utils/StringUtils.h>
 #include <itkImageFileReader.h>
 #include <itkOrientImageFilter.h>
@@ -912,7 +913,7 @@ void Shape::set_reconstruction_transforms(std::vector<vtkSmartPointer<vtkTransfo
 
 //---------------------------------------------------------------------------
 void Shape::load_feature_from_scalar_file(std::string filename, std::string feature_name) {
-  if (!boost::filesystem::exists(filename)) {
+  if (!ShapeWorksUtils::file_exists(filename)) {
     return;
   }
 

Libs/Common/ShapeworksUtils.cpp

@@ -1,11 +1,19 @@
+// boost::nowide's utf8_codecvt is header-only, but including it pulls in the auto-link pragma
+// that asks MSVC for boost_nowide.lib, which the dependency build does not produce.  This must
+// precede every Boost header.
+#define BOOST_NOWIDE_NO_LIB
+
 #include "ShapeworksUtils.h"
 
 #include <Logging.h>
 #include <sys/stat.h>
 #include <tbb/global_control.h>
 #include <tbb/info.h>
+#include <vtksys/SystemTools.hxx>
 
 #include <boost/filesystem.hpp>
+#include <boost/nowide/utf8_codecvt.hpp>
+#include <locale>
 
 namespace shapeworks {
 
@@ -19,16 +27,33 @@ void ShapeWorksUtils::set_rng_seed(const unsigned seed) {
   mt_.seed(rng_seed_);
 }
 
+//-----------------------------------------------------------------------------
+void ShapeWorksUtils::initialize_path_handling() {
+  // boost::filesystem stores a path in the character type the OS API uses: char on POSIX, but
+  // wchar_t on Windows.  There every narrow string is converted using the codecvt facet imbued
+  // into boost::filesystem::path, which by default comes from the global locale rather than being
+  // UTF-8.  The rest of ShapeWorks (Qt, VTK, ITK) speaks UTF-8, so without this a path holding
+  // non-ASCII characters is mangled on the way in and no longer resolves.  On POSIX no conversion
+  // takes place and this only affects the wide-string accessors.
+  boost::filesystem::path::imbue(std::locale(std::locale(), new boost::nowide::utf8_codecvt<wchar_t>));
+}
+
 //-----------------------------------------------------------------------------
 bool ShapeWorksUtils::is_directory(const std::string& pathname) {
-  boost::system::error_code ec;
-  return boost::filesystem::is_directory(pathname, ec);
+  return vtksys::SystemTools::FileIsDirectory(pathname);
 }
 
 //-----------------------------------------------------------------------------
 bool ShapeWorksUtils::file_exists(const std::string& filename) {
-  boost::system::error_code ec;
-  return boost::filesystem::is_regular_file(filename, ec);
+  // vtksys rather than boost::filesystem: on Windows it decodes UTF-8 and applies the \\?\ prefix,
+  // so it agrees with the VTK and ITK readers about which paths exist.  boost::filesystem does
+  // neither, and rejects long or non-ASCII paths that those readers open without complaint.
+  return vtksys::SystemTools::FileExists(filename, true);
+}
+
+//-----------------------------------------------------------------------------
+bool ShapeWorksUtils::path_exists(const std::string& pathname) {
+  return vtksys::SystemTools::PathExists(pathname);
 }
 
 //-----------------------------------------------------------------------------

Libs/Common/ShapeworksUtils.h

@@ -31,14 +31,19 @@ class ShapeWorksUtils {
   /// generates random number
   static unsigned get_random_number() { return mt_(); }
 
+  //! Configure process-wide filesystem path handling.  Call this once at startup, before any path
+  //! is constructed from a narrow string.
+  static void initialize_path_handling();
+
   /// returns true if pathname is a directory
-  // TODO: in C++17 this is a standard function
   static bool is_directory(const std::string& pathname);
 
-  /// returns true if filename exists
-  // TODO: in C++17 this is a standard function
+  /// returns true if filename exists and can be read
   static bool file_exists(const std::string& filename);
 
+  /// returns true if a file or a directory exists at pathname
+  static bool path_exists(const std::string& pathname);
+
   //! Set up the console logging options
   static void setup_console_logging(bool show_progress, bool xml_status);
 

Libs/Image/Image.cpp

@@ -117,7 +117,7 @@ Image::ImageType::Pointer Image::read(const std::string& pathname) {
   }
 
   // check if it exists
-  if (!boost::filesystem::exists(pathname)) {
+  if (!ShapeWorksUtils::file_exists(pathname)) {
     throw std::runtime_error("File does not exist: " + pathname);
   }
 

Libs/Mesh/Mesh.cpp

@@ -5,7 +5,9 @@
 #include <tbb/parallel_for.h>
 
 #include <algorithm>
+#include <numeric>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 // libigl
@@ -20,6 +22,7 @@
 #include <vtkAlgorithm.h>
 #include <vtkAppendPolyData.h>
 #include <vtkButterflySubdivisionFilter.h>
+#include <vtkCellArrayIterator.h>
 #include <vtkCenterOfMass.h>
 #include <vtkCleanPolyData.h>
 #include <vtkClipClosedSurface.h>
@@ -29,6 +32,7 @@
 #include <vtkFeatureEdges.h>
 #include <vtkFillHolesFilter.h>
 #include <vtkGenericCell.h>
+#include <vtkGeometryFilter.h>
 #include <vtkGradientFilter.h>
 #include <vtkImageData.h>
 #include <vtkImageStencil.h>
@@ -57,6 +61,7 @@
 #include <vtkSmoothPolyDataFilter.h>
 #include <vtkStaticCellLocator.h>
 #include <vtkTransformPolyDataFilter.h>
+#include <vtkUnstructuredGridReader.h>
 #include <vtkWindowedSincPolyDataFilter.h>
 #include <vtkXMLPolyDataReader.h>
 #include <vtkXMLPolyDataWriter.h>
@@ -113,6 +118,78 @@ vtkSmartPointer<vtkPolyData> update_and_check(vtkAlgorithm* reader, const std::s
   }
   return vtkPolyData::SafeDownCast(reader->GetOutputDataObject(0));
 }
+
+//! Remove polygons that appear more than once (same point ids, in any order or winding). Duplicated
+//! faces typically form zero-thickness fins hanging off the surface, which create non-manifold edges
+//! and corrupt normals. All copies are removed unless that would open a hole in the surface, in which
+//! case one copy is kept. Returns the number of faces removed.
+vtkIdType remove_duplicate_faces(vtkPolyData* poly_data) {
+  auto polys = poly_data->GetPolys();
+  const vtkIdType num_polys = polys->GetNumberOfCells();
+  if (num_polys == 0) {
+    return 0;
+  }
+
+  auto edge_key = [](vtkIdType a, vtkIdType b) {
+    return (static_cast<uint64_t>(std::min(a, b)) << 32) | static_cast<uint64_t>(std::max(a, b));
+  };
+
+  std::vector<std::vector<vtkIdType>> faces(num_polys);
+  std::unordered_map<uint64_t, int> edge_count;
+  auto iter = vtk::TakeSmartPointer(polys->NewIterator());
+  vtkIdType index = 0;
+  for (iter->GoToFirstCell(); !iter->IsDoneWithTraversal(); iter->GoToNextCell(), index++) {
+    vtkIdType num_points;
+    const vtkIdType* points;
+    iter->GetCurrentCell(num_points, points);
+    faces[index].assign(points, points + num_points);
+    for (vtkIdType i = 0; i < num_points; i++) {
+      edge_count[edge_key(points[i], points[(i + 1) % num_points])]++;
+    }
+    std::sort(faces[index].begin(), faces[index].end());
+  }
+
+  std::vector<vtkIdType> order(num_polys);
+  std::iota(order.begin(), order.end(), 0);
+  std::sort(order.begin(), order.end(), [&](vtkIdType a, vtkIdType b) { return faces[a] < faces[b]; });
+
+  // polys follow verts and lines in the cell id ordering
+  const vtkIdType poly_offset = poly_data->GetNumberOfVerts() + poly_data->GetNumberOfLines();
+  std::vector<vtkIdType> to_delete;
+  for (vtkIdType start = 0; start < num_polys;) {
+    vtkIdType end = start + 1;
+    while (end < num_polys && faces[order[end]] == faces[order[start]]) {
+      end++;
+    }
+    const int copies = end - start;
+    if (copies > 1) {
+      // if removing every copy would leave an edge with a single face, this duplicates a real surface face
+      const auto& face = faces[order[start]];
+      bool opens_hole = false;
+      for (size_t i = 0; i < face.size(); i++) {
+        for (size_t j = i + 1; j < face.size(); j++) {
+          auto it = edge_count.find(edge_key(face[i], face[j]));
+          if (it != edge_count.end() && it->second - copies == 1) {
+            opens_hole = true;
+          }
+        }
+      }
+      for (vtkIdType k = opens_hole ? start + 1 : start; k < end; k++) {
+        to_delete.push_back(poly_offset + order[k]);
+      }
+    }
+    start = end;
+  }
+
+  if (!to_delete.empty()) {
+    poly_data->BuildCells();
+    for (auto id : to_delete) {
+      poly_data->DeleteCell(id);
+    }
+    poly_data->RemoveDeletedCells();
+  }
+  return to_delete.size();
+}
 }  // namespace
 
 vtkSmartPointer<vtkPolyData> MeshReader::read(const std::string& pathname) {
@@ -125,6 +202,17 @@ vtkSmartPointer<vtkPolyData> MeshReader::read(const std::string& pathname) {
 
   try {
     if (StringUtils::hasSuffix(pathname, ".vtk")) {
+      // legacy VTK files may store a surface as an unstructured grid; extract its surface as polydata
+      auto ug_reader = vtkSmartPointer<vtkUnstructuredGridReader>::New();
+      ug_reader->SetFileName(pathname.c_str());
+      if (ug_reader->IsFileUnstructuredGrid()) {
+        ug_reader->SetReadAllScalars(1);
+        update_and_check(ug_reader, pathname);
+        auto geometry = vtkSmartPointer<vtkGeometryFilter>::New();
+        geometry->SetInputConnection(ug_reader->GetOutputPort());
+        return update_and_check(geometry, pathname);
+      }
+
       auto reader = vtkSmartPointer<vtkPolyDataReader>::New();
       reader->SetFileName(pathname.c_str());
       reader->SetReadAllScalars(1);
@@ -390,7 +478,24 @@ Mesh& Mesh::remesh(int numVertices, double adaptivity) {
 
   numVertices = std::max<int>(numVertices, 1);
 
-  remesh->SetForceManifold(true);
+  // ACVD's ForceManifold treats clusters on an open boundary as non-manifold and keeps splitting
+  // them, piling extra vertices and sliver triangles along the boundary (e.g. a mesh clipped by a
+  // cutting plane ends up with more vertices than requested). Only force manifold on closed meshes.
+  // geometry only: vtkFeatureEdges copies point data, which reads out of bounds on a field that
+  // doesn't have one value per point
+  auto geometry = vtkSmartPointer<vtkPolyData>::New();
+  geometry->SetPoints(this->poly_data_->GetPoints());
+  geometry->SetPolys(this->poly_data_->GetPolys());
+  auto boundary_edges = vtkSmartPointer<vtkFeatureEdges>::New();
+  boundary_edges->SetInputData(geometry);
+  boundary_edges->BoundaryEdgesOn();
+  boundary_edges->FeatureEdgesOff();
+  boundary_edges->NonManifoldEdgesOff();
+  boundary_edges->ManifoldEdgesOff();
+  boundary_edges->Update();
+  bool is_closed = boundary_edges->GetOutput()->GetNumberOfCells() == 0;
+
+  remesh->SetForceManifold(is_closed);
   remesh->SetInput(surf);
   remesh->SetFileLoadSaveOption(0);
   remesh->SetConsoleOutput(0);
@@ -499,14 +604,22 @@ Mesh& Mesh::fillHoles(double hole_size) {
 
 Mesh& Mesh::clean() {
   TIME_SCOPE("Mesh::clean");
-  auto clean = vtkSmartPointer<vtkCleanPolyData>::New();
-  clean->ConvertPolysToLinesOff();
-  clean->ConvertLinesToPointsOff();
-  clean->ConvertStripsToPolysOff();
-  clean->PointMergingOn();
-  clean->SetInputData(poly_data_);
-  clean->Update();
-  poly_data_ = clean->GetOutput();
+  auto run_clean = [](vtkPolyData* input) {
+    auto clean = vtkSmartPointer<vtkCleanPolyData>::New();
+    clean->ConvertPolysToLinesOff();
+    clean->ConvertLinesToPointsOff();
+    clean->ConvertStripsToPolysOff();
+    clean->PointMergingOn();
+    clean->SetInputData(input);
+    clean->Update();
+    return vtkSmartPointer<vtkPolyData>(clean->GetOutput());
+  };
+
+  poly_data_ = run_clean(poly_data_);
+  if (remove_duplicate_faces(poly_data_) > 0) {
+    // drop points left unused by the removed faces
+    poly_data_ = run_clean(poly_data_);
+  }
   invalidateLocators();
   return *this;
 }

Libs/Mesh/Mesh.h

@@ -104,7 +104,7 @@ class Mesh {
   /// finds holes in a mesh and closes them
   Mesh& fillHoles(double hole_size = 1000.0);
 
-  /// clean mesh
+  /// clean mesh: merges duplicate points, removes unused points, and removes duplicated faces
   Mesh& clean();
 
   /// samples image data values at point locations specified by image