Propagation

.gitignore

@@ -7,3 +7,5 @@ CMakeLists.txt.user.*
 .8.un~
 docs/_build
 *.asv
+build/
+.vscode/

CMakeLists.txt

@@ -67,7 +67,17 @@ FIND_PACKAGE(ITK 5.2.0 REQUIRED)
 INCLUDE(${ITK_USE_FILE})
 
 # VTK - required for quality mesh transformations
-FIND_PACKAGE(VTK 9.1.0 REQUIRED COMPONENTS CommonCore IOCore IOLegacy IOPLY IOGeometry FiltersModeling)
+FIND_PACKAGE(VTK 9.1.0 REQUIRED COMPONENTS
+  CommonCore
+  IOCore
+  IOLegacy
+  IOPLY
+  IOGeometry
+  IOImage
+  IOXML
+  FiltersCore
+  FiltersGeneral 
+  FiltersModeling)
 SET(GREEDY_VTK_LIBRARIES ${VTK_LIBRARIES})
 
 # Deal with FFTW - only used by experimental LDDMM code
@@ -240,8 +250,15 @@ SET(GREEDY_API_LIBS greedyapi
   ${FFTWF_LIB} ${FFTWF_THREADS_LIB}
   ${SPARSE_LIBRARY})
 
+# propagation api
+add_subdirectory(src/propagation propagation)
+
+SET(PROPAGATION_SRC
+  src/propagation/greedy_propagation.cxx
+)
+
 # List of installable targets
-SET(GREEDY_INSTALL_BIN_TARGETS greedy greedy_template_average multi_chunk_greedy)
+SET(GREEDY_INSTALL_BIN_TARGETS greedy greedy_template_average multi_chunk_greedy greedy_propagation)
 SET(GREEDY_INSTALL_LIB_TARGETS greedyapi)
 
 
@@ -300,8 +317,13 @@ IF(BUILD_CLI)
   ADD_EXECUTABLE(test_greedy testing/src/GreedyTestDriver.cxx)
   TARGET_LINK_LIBRARIES(test_greedy ${GREEDY_API_LIBS})
 
+  ADD_EXECUTABLE(test_propagation testing/src/propagation/propagation_test.cxx)
+  TARGET_LINK_LIBRARIES(test_propagation propagationapi)
+  TARGET_INCLUDE_DIRECTORIES(test_propagation PUBLIC ${GREEDY_SOURCE_DIR}/src/propagation )
+
   ADD_EXECUTABLE(multi_chunk_greedy ${CHUNK_GREEDY_SRC})
   TARGET_LINK_LIBRARIES(multi_chunk_greedy ${GREEDY_API_LIBS})
+
 ENDIF(BUILD_CLI)
 
 # Install command-line executables
@@ -421,6 +443,9 @@ IF(NOT GREEDY_BUILD_AS_SUBPROJECT)
   ADD_TEST(NAME "Phantom_NCC_Sim_NoMask"  COMMAND test_greedy phantom 1 3 NCC 7 0 WORKING_DIRECTORY ${TESTING_DATADIR})
   ADD_TEST(NAME "Phantom_WNCC_Sim_NoMask" COMMAND test_greedy phantom 1 3 WNCC 7 0 WORKING_DIRECTORY ${TESTING_DATADIR})
 
+  # Tests for greedy_propagation
+  ADD_TEST(NAME "propagation_basic" COMMAND test_propagation basic WORKING_DIRECTORY ${TESTING_DATADIR})
+  ADD_TEST(NAME "propagation_extra_mesh" COMMAND test_propagation extra_mesh WORKING_DIRECTORY ${TESTING_DATADIR})
   # Add tests for lmshoot
   IF(GREEDY_BUILD_LMSHOOT)
   ADD_TEST(NAME "lmshoot_regression" COMMAND lmshoot_test shoot_regression.mat WORKING_DIRECTORY ${TESTING_DATADIR}/lmshoot)

src/CommandLineHelper.h

@@ -256,6 +256,21 @@ class CommandLineHelper
     return val;
   }
 
+  /**
+   * Read output directory path
+   */
+  std::string read_output_dir()
+  {
+    std::string dir = read_arg();
+    if(this->data_root.length())
+      dir = itksys::SystemTools::CollapseFullPath(dir, data_root);
+
+    if(!itksys::SystemTools::PathExists(dir.c_str()))
+      throw GreedyException("Folder '%s' does not exist", dir.c_str());
+
+    return dir;
+  }
+
   /**
    * Check if a string ends with another string and return the
    * substring without the suffix
@@ -360,27 +375,26 @@ class CommandLineHelper
     return vector;
   }
 
-  std::vector<int> read_int_vector()
+  std::vector<int> read_int_vector(char delimiter = 'x')
   {
     std::string arg = read_arg();
     std::istringstream f(arg);
     std::string s;
     std::vector<int> vector;
-    while (getline(f, s, 'x'))
+    while (getline(f, s, delimiter))
       {
       errno = 0; char *pend;
       long val = std::strtol(s.c_str(), &pend, 10);
 
       if(errno || *pend)
-        throw GreedyException("Expected an integer vector as parameter to '%s', instead got '%s'",
-                              current_command.c_str(), arg.c_str());
+        throw GreedyException("Expected an integer vector delimited by '%c' as parameter to '%s', instead got '%s'",
+                              delimiter, current_command.c_str(), arg.c_str());
       vector.push_back((int) val);
       }
 
     if(!vector.size())
-      throw GreedyException("Expected an integer vector as parameter to '%s', instead got '%s'",
-                            current_command.c_str(), arg.c_str());
-
+      throw GreedyException("Expected an integer vector delimited by '%c' as parameter to '%s', instead got '%s'",
+                            delimiter, current_command.c_str(), arg.c_str());
     return vector;
   }
 

src/GreedyAPI.cxx

@@ -416,6 +416,28 @@ ::ReadImageViaCache(const std::string &filename,
   return pointer;
 }
 
+template <unsigned int VDim, typename TReal>
+typename GreedyApproach<VDim, TReal>::MeshPointer
+GreedyApproach<VDim, TReal>
+::ReadMeshViaCache(const std::string &filename)
+{
+  typename MeshCache::const_iterator it = m_MeshCache.find(filename);
+  if(it != m_MeshCache.cend())
+    {
+    vtkObject *cached_object = it->second.target;
+    MeshType *mesh = dynamic_cast<MeshType*>(cached_object);
+    if (!mesh)
+      throw GreedyException("Cached mesh %s cannot be cast to type %s",
+                            filename.c_str(), typeid(MeshType).name());
+    MeshPointer pMesh = DeepCopyMesh(mesh); // important to avoid in-place mutation
+    return pMesh;
+    }
+
+  // Read the mesh using mesh io reader
+  return ReadMesh(filename.c_str());
+}
+
+
 template <unsigned int VDim, typename TReal>
 template <class TObject>
 TObject *
@@ -542,6 +564,27 @@ ::WriteImageViaCache(TImage *img, const std::string &filename, itk::IOComponentE
     }
 }
 
+template <unsigned int VDim, typename TReal>
+void
+GreedyApproach<VDim, TReal>
+::WriteMeshViaCache(MeshType *mesh, const std::string &filename)
+{
+  typename MeshCache::const_iterator it = m_MeshCache.find(filename);
+  if (it != m_MeshCache.end())
+    {
+    auto *cached = dynamic_cast<MeshType*>(it->second.target);
+    if (!cached)
+      throw GreedyException("Cached mesh %s cannot be cast to %s",
+                            filename.c_str(), typeid(MeshType*).name());
+    cached->DeepCopy(mesh);
+    }
+
+  if (it == m_MeshCache.end() || it->second.force_write)
+    {
+    WriteMesh(mesh, filename.c_str());
+    }
+}
+
 
 
 #include <itkBinaryErodeImageFilter.h>
@@ -1719,7 +1762,7 @@ ::RunDeformable(GreedyParameters &param)
   if(param.tjr_param.weight > 0.0)
     {
     // Read the mesh
-    vtkSmartPointer<vtkPointSet> point_set = ReadMesh(param.tjr_param.tetra_mesh.c_str());
+    vtkSmartPointer<vtkPointSet> point_set = ReadMeshViaCache(param.tjr_param.tetra_mesh.c_str());
     vtkSmartPointer<vtkUnstructuredGrid> tetra = dynamic_cast<vtkUnstructuredGrid *>(point_set.GetPointer());
     if(!tetra)
       throw GreedyException("Mesh %s is not an UnstructuredGrid!", param.tjr_param.tetra_mesh.c_str());
@@ -1905,8 +1948,9 @@ ::RunDeformable(GreedyParameters &param)
       tm_Gradient.Stop();
 
       // Print a report for this iteration
-      std::cout << this->PrintIter(level, iter, metric_report, reg_report) << std::endl;
-      fflush(stdout);
+      std::string iter_line = this->PrintIter(level, iter, metric_report, reg_report);
+      gout.printf("%s\n", iter_line.c_str());
+      gout.flush();
 
       // Record the metric value in the log
       this->RecordMetricValue(metric_report);
@@ -2279,7 +2323,7 @@ ::RunDeformableOptimization(GreedyParameters &param)
   if(param.tjr_param.weight > 0.0)
     {
     // Read the mesh
-    vtkSmartPointer<vtkPointSet> point_set = ReadMesh(param.tjr_param.tetra_mesh.c_str());
+    vtkSmartPointer<vtkPointSet> point_set = ReadMeshViaCache(param.tjr_param.tetra_mesh.c_str());
     vtkSmartPointer<vtkUnstructuredGrid> tetra = dynamic_cast<vtkUnstructuredGrid *>(point_set.GetPointer());
     if(!tetra)
       throw GreedyException("Mesh %s is not an UnstructuredGrid!", param.tjr_param.tetra_mesh.c_str());
@@ -2451,8 +2495,9 @@ ::RunDeformableOptimization(GreedyParameters &param)
         }
 
       // Print a report for this iteration
-      std::cout << this->PrintIter(level, iter, metric_report, reg_report) << std::endl;
-      fflush(stdout);
+      std::string iter_line = this->PrintIter(level, iter, metric_report, reg_report);
+      gout.printf("%s\n", iter_line.c_str());
+      gout.flush();
 
       // Record the metric value in the log
       this->RecordMetricValue(metric_report);
@@ -3251,7 +3296,7 @@ ::RunReslice(GreedyParameters &param)
   std::vector<MeshPointer> meshes, original_meshes;
   for(unsigned int i = 0; i < r_param.meshes.size(); i++)
     {
-    vtkSmartPointer<vtkPointSet> mesh = ReadMesh(r_param.meshes[i].fixed.c_str());
+    vtkSmartPointer<vtkPointSet> mesh = ReadMeshViaCache(r_param.meshes[i].fixed.c_str());
     meshes.push_back(mesh);
 
     if(r_param.meshes[i].jacobian_mode)
@@ -3425,7 +3470,7 @@ ::RunReslice(GreedyParameters &param)
     if(r_param.meshes[i].jacobian_mode)
       WriteJacobianMesh(original_meshes[i], meshes[i], r_param.meshes[i].output.c_str());
     else
-      WriteMesh(meshes[i], r_param.meshes[i].output.c_str());
+      WriteMeshViaCache(meshes[i], r_param.meshes[i].output.c_str());
     }
 
 
@@ -3739,6 +3784,15 @@ ::AddCachedInputObject(std::string key, itk::Object *object)
   m_ImageCache[key].force_write = false;
 }
 
+template <unsigned int VDim, typename TReal>
+void GreedyApproach<VDim, TReal>
+::AddCachedInputObject(std::string key, vtkObject *object)
+{
+  m_MeshCache[key].target = object;
+  m_MeshCache[key].force_write = false;
+}
+
+
 template <unsigned int VDim, typename TReal>
 void GreedyApproach<VDim, TReal>
 ::AddCachedOutputObject(std::string key, itk::Object *object, bool force_write)
@@ -3768,6 +3822,14 @@ ::GetCachedObjectNames() const
   return keys;
 }
 
+template <unsigned int VDim, typename TReal>
+void GreedyApproach<VDim, TReal>
+::AddCachedOutputObject(std::string key, vtkObject *object, bool force_write)
+{
+  m_MeshCache[key].target = object;
+  m_MeshCache[key].force_write = force_write;
+}
+
 template <unsigned int VDim, typename TReal>
 const typename GreedyApproach<VDim,TReal>::MetricLogType &
 GreedyApproach<VDim,TReal>

src/GreedyAPI.h

@@ -37,6 +37,7 @@
 #include <map>
 #include "itkCommand.h"
 #include <vtkSmartPointer.h>
+#include <vtkObject.h>
 
 template <typename T, unsigned int V> class MultiImageOpticalFlowHelper;
 
@@ -85,7 +86,8 @@ class GreedyApproach
   };
 
   // Mesh data structures
-  typedef vtkSmartPointer<vtkPointSet> MeshPointer;
+  typedef vtkPointSet MeshType;
+  typedef vtkSmartPointer<MeshType> MeshPointer;
   typedef std::vector<MeshPointer> MeshArray;
 
   static void ConfigThreads(const GreedyParameters &param);
@@ -141,6 +143,7 @@ class GreedyApproach
    *
    */
   void AddCachedInputObject(std::string key, itk::Object *object);
+  void AddCachedInputObject(std::string key, vtkObject *object);
 
   /**
    * Add an image/matrix to the output cache. This has the same behavior as
@@ -153,6 +156,7 @@ class GreedyApproach
    * will be allocated. It can then me accessed using GetCachedObject()
    */
   void AddCachedOutputObject(std::string key, itk::Object *object, bool force_write = false);
+  void AddCachedOutputObject(std::string key, vtkObject *object, bool force_write = false);
 
   /**
    * Get a cached object by name
@@ -291,9 +295,17 @@ class GreedyApproach
     bool force_write;
   };
 
+  struct VTKCacheEntry {
+    vtkObject *target;
+    bool force_write;
+  };
+
   typedef std::map<std::string, CacheEntry> ImageCache;
   ImageCache m_ImageCache;
 
+  typedef std::map<std::string, VTKCacheEntry> MeshCache;
+  MeshCache m_MeshCache;
+
   // A log of metric values used during registration - so metric can be looked up
   // in the callbacks to RunAffine, etc.
   MetricLogType m_MetricLog;
@@ -311,6 +323,8 @@ class GreedyApproach
   itk::SmartPointer<TImage> ReadImageViaCache(const std::string &filename,
                                               itk::IOComponentEnum *comp_type = NULL);
 
+  MeshPointer ReadMeshViaCache(const std::string &filename);
+
   template<class TObject> TObject *CheckCache(const std::string &filename) const;
 
   // Get a filename for dumping intermediate outputs
@@ -330,6 +344,8 @@ class GreedyApproach
   void WriteCompressedWarpInPhysicalSpaceViaCache(
     ImageBaseType *moving_ref_space, VectorImageType *warp, const char *filename, double precision);
 
+  void WriteMeshViaCache(MeshType *mesh, const std::string &filename);
+
   // Compute the moments of a composite image (mean and covariance matrix of coordinate weighted by intensity)
   void ComputeImageMoments(CompositeImageType *image, const vnl_vector<float> &weights, VecFx &m1, MatFx &m2);
 

src/GreedyMeshIO.cxx

@@ -12,6 +12,8 @@
 #include <vtkBYUReader.h>
 #include <vtkBYUWriter.h>
 #include <vtkOBJReader.h>
+#include <vtkXMLPolyDataReader.h>
+#include <vtkXMLPolyDataWriter.h>
 #include <vtkTetra.h>
 #include <vtkTriangle.h>
 #include <vtkDoubleArray.h>
@@ -73,6 +75,8 @@ vtkSmartPointer<TMesh> ReadMeshByExtension(const char *fname)
     else
       throw GreedyException("No mesh reader for file %s", fname);
     }
+    else if(fn_str.rfind(".vtp") == fn_str.length() - 4)
+      return ReadMesh<vtkXMLPolyDataReader, TMesh>(fname);
   else
     throw GreedyException("No mesh reader for file %s", fname);
 }
@@ -96,6 +100,11 @@ void WriteMeshByExtension(TMesh *mesh, const char *fname)
     else if (usg)
       WriteMesh<vtkUnstructuredGridWriter, vtkUnstructuredGrid>(usg, fname);
     }
+  else if (fn_str.rfind(".vtp") == fn_str.length() - 4)
+    {
+    vtkPolyData *pd = dynamic_cast<vtkPolyData *>(mesh);
+    WriteMesh<vtkXMLPolyDataWriter, TMesh>(pd, fname);
+    }
   else
     throw GreedyException("No mesh writer for file %s", fname);
 }