Merge v6.6.1 into master

Libs/Analyze/MeshGenerator.cpp

@@ -7,6 +7,7 @@
 #include <QMeshWarper.h>
 #include <SurfaceReconstructor.h>
 #include <Utils/StringUtils.h>
+#include <itkConstantPadImageFilter.h>
 #include <itkImageFileReader.h>
 #include <itkOrientImageFilter.h>
 #include <itkPoint.h>
@@ -16,7 +17,6 @@
 #include <vtkPolyDataNormals.h>
 
 #include <QFileInfo>
-#include <limits>
 
 namespace shapeworks {
 
@@ -99,12 +99,16 @@ MeshHandle MeshGenerator::build_mesh_from_image(ImageType::Pointer image, float
 
   try {
     // only interested in 1's and 0's
-    Image itk_image = Image(image);
-    if (!itk_image.isDistanceTransform()) {
-      itk_image.binarize(0, 1);
-      image = itk_image.getITKImage();
+    Image sw_image = Image(image);
+    if (!sw_image.isDistanceTransform()) {
+      sw_image.binarize();
+      image = sw_image.getITKImage();
     }
 
+    // pad the image in case the segmentation is on the edge
+    sw_image.pad(1);
+    image = sw_image.getITKImage();
+
     // connect to VTK
     vtkSmartPointer<vtkImageImport> vtk_image = vtkSmartPointer<vtkImageImport>::New();
     itk::VTKImageExport<ImageType>::Pointer itk_exporter = itk::VTKImageExport<ImageType>::New();
@@ -165,6 +169,8 @@ MeshHandle MeshGenerator::build_mesh_from_file(std::string filename, float iso_v
     }
   } else if (is_image) {
     try {
+      ImageUtils::register_itk_factories();
+
       // read file using ITK
       using ReaderType = itk::ImageFileReader<ImageType>;
       ReaderType::Pointer reader = ReaderType::New();

Libs/Analyze/MeshWorker.cpp

@@ -19,6 +19,9 @@ void MeshWorker::run()
 {
   // build the mesh using our MeshGenerator
   auto item = this->queue_->get_next_work_item();
+  if (!item) {
+    return;
+  }
 
   MeshHandle mesh = this->mesh_generator_->build_mesh(*item);
 

Libs/Analyze/Shape.cpp

@@ -125,10 +125,8 @@ void Shape::recompute_original_surface() {
   original_meshes_.set_mesh(0, mesh_handle);
 }
 
-
 //---------------------------------------------------------------------------
-void Shape::ensure_segmentation()
-{
+void Shape::ensure_segmentation() {
   if (get_segmentation()) {
     return;
   }
@@ -159,8 +157,7 @@ void Shape::ensure_segmentation()
 //---------------------------------------------------------------------------
 MeshGroup Shape::get_groomed_meshes(bool wait) {
   if (!subject_) {
-    std::cerr << "Error: asked for groomed meshes when none are present!\n";
-    assert(0);
+    return {};
   }
 
   if (!groomed_meshes_.valid() || groomed_meshes_.meshes().size() != subject_->get_number_of_domains()) {

Libs/Common/Logging.cpp

@@ -13,11 +13,11 @@ namespace shapeworks {
 static std::string create_header(const int line, const char* filename, const char* function = "") {
   const char* name = (strrchr(filename, '/') ? strrchr(filename, '/') + 1 : filename);
   const char* name2 = (strrchr(name, '\\') ? strrchr(name, '\\') + 1 : name);
-  const char* function_name = (strrchr(function, ':') ? strrchr(function, ':') + 1 : function);
   if (!function) {
     std::string header = "[" + std::string(name2) + "|" + std::to_string(line) + "]";
     return header;
   } else {
+    const char* function_name = (strrchr(function, ':') ? strrchr(function, ':') + 1 : function);
     std::string header = "[" + std::string(name2) + "|" + std::string(function_name) + "|" + std::to_string(line) + "]";
     return header;
   }
@@ -69,6 +69,14 @@ void Logging::log_message(const std::string& message, const int line, const char
   }
 }
 
+//-----------------------------------------------------------------------------
+void Logging::log_only(const std::string& message, const int line, const char* file) const {
+  spd::info(message);
+  if (log_open_) {
+    spd::get("file")->info(message);
+  }
+}
+
 //-----------------------------------------------------------------------------
 void Logging::log_stack(const std::string& message) const {
   spd::error(message);

Libs/Common/Logging.h

@@ -102,6 +102,9 @@ class Logging {
   //! Log a message, use SW_LOG macro
   void log_message(const std::string& message, const int line, const char* file) const;
 
+  //! Log a message, use SW_LOG_ONLY macro
+  void log_only(const std::string& message, const int line, const char* file) const;
+
   //! Log a stack trace message, use SW_LOG_STACK macro
   void log_stack(const std::string& message) const;
 
@@ -171,6 +174,10 @@ class Logging {
 #define SW_LOG(message, ...) \
   shapeworks::Logging::Instance().log_message(safe_format(message, ##__VA_ARGS__), __LINE__, __FILE__);
 
+//! Log only macro
+#define SW_LOG_ONLY(message, ...) \
+  shapeworks::Logging::Instance().log_only(safe_format(message, ##__VA_ARGS__), __LINE__, __FILE__);
+
 //! Log warning macro
 #define SW_WARN(message, ...) \
   shapeworks::Logging::Instance().log_warning(safe_format(message, ##__VA_ARGS__), __LINE__, __FILE__)

Libs/Groom/GroomParameters.h

@@ -137,7 +137,6 @@ class GroomParameters {
   bool get_skip_grooming();
   void set_skip_grooming(bool skip);
 
-
   bool get_shared_boundary();
   void set_shared_boundary(bool shared_boundary);
 
@@ -152,6 +151,8 @@ class GroomParameters {
 
   void restore_defaults();
 
+  Parameters get_parameters() const { return params_; }
+
   // constants
   const static std::string GROOM_SMOOTH_VTK_LAPLACIAN_C;
   const static std::string GROOM_SMOOTH_VTK_WINDOWED_SINC_C;

Libs/Image/Image.cpp

@@ -44,26 +44,8 @@
 #include "ShapeworksUtils.h"
 #include "itkTPGACLevelSetImageFilter.h"  // actually a shapeworks class, not itk
 
-// ITK image factories
-#include <itkMetaImageIOFactory.h>
-#include <itkNiftiImageIOFactory.h>
-#include <itkNrrdImageIOFactory.h>
-
 namespace shapeworks {
 
-namespace {
-void register_factories() {
-  static bool registered = false;
-  if (!registered) {
-    // register all the factories
-    itk::NrrdImageIOFactory::RegisterOneFactory();
-    itk::NiftiImageIOFactory::RegisterOneFactory();
-    itk::MetaImageIOFactory::RegisterOneFactory();
-    registered = true;
-  }
-}
-}  // namespace
-
 Image::Image(const Dims dims) : itk_image_(ImageType::New()) {
   ImageType::RegionType region;
   region.SetSize(dims);
@@ -119,7 +101,7 @@ Image& Image::operator=(Image&& img) {
 }
 
 Image::ImageType::Pointer Image::read(const std::string& pathname) {
-  register_factories();
+  ImageUtils::register_itk_factories();
 
   if (pathname.empty()) {
     throw std::invalid_argument("Empty pathname");

Libs/Image/ImageUtils.cpp

@@ -2,6 +2,11 @@
 
 #include <itkPointSet.h>
 
+// ITK image factories
+#include <itkMetaImageIOFactory.h>
+#include <itkNiftiImageIOFactory.h>
+#include <itkNrrdImageIOFactory.h>
+
 namespace shapeworks {
 
 PhysicalRegion ImageUtils::boundingBox(const std::vector<std::string>& filenames, Image::PixelType isoValue) {
@@ -39,6 +44,18 @@ PhysicalRegion ImageUtils::boundingBox(const std::vector<std::reference_wrapper<
   return bbox;
 }
 
+void ImageUtils::register_itk_factories()
+{
+  static bool registered = false;
+  if (!registered) {
+    // register all the factories
+    itk::NrrdImageIOFactory::RegisterOneFactory();
+    itk::NiftiImageIOFactory::RegisterOneFactory();
+    itk::MetaImageIOFactory::RegisterOneFactory();
+    registered = true;
+  }
+}
+
 ImageUtils::TPSTransform::Pointer ImageUtils::createWarpTransform(const std::string& source_landmarks_file,
                                                                   const std::string& target_landmarks_file,
                                                                   const int stride) {

Libs/Image/ImageUtils.h

@@ -22,6 +22,9 @@ class ImageUtils {
   using TPSTransform = itk::ThinPlateSplineKernelTransform<double, 3>;
   static TPSTransform::Pointer createWarpTransform(const std::string& source_landmarks_file,
                                                    const std::string& target_landmarks_file, const int stride = 1);
+
+  static void register_itk_factories();
+
 };
 
 }  // namespace shapeworks

Libs/Mesh/MeshWarper.cpp

@@ -56,6 +56,9 @@ vtkSmartPointer<vtkPolyData> MeshWarper::build_mesh(const Eigen::MatrixXd& parti
 //---------------------------------------------------------------------------
 void MeshWarper::set_reference_mesh(vtkSmartPointer<vtkPolyData> reference_mesh,
                                     const Eigen::MatrixXd& reference_particles, const Eigen::MatrixXd& landmarks) {
+  // lock so that we don't swap out the reference mesh while we are using it
+  std::scoped_lock lock(mutex);
+
   if (this->incoming_reference_mesh_ == reference_mesh) {
     if (this->reference_particles_.size() == reference_particles.size()) {
       if (this->reference_particles_ == reference_particles && landmarks_points_ == landmarks) {

Libs/Optimize/OptimizeParameters.cpp

@@ -864,3 +864,6 @@ double OptimizeParameters::get_shared_boundary_weight() { return params_.get(Key
 
 //---------------------------------------------------------------------------
 void OptimizeParameters::set_shared_boundary_weight(double value) { params_.set(Keys::shared_boundary_weight, value); }
+
+//---------------------------------------------------------------------------
+Parameters OptimizeParameters::get_parameters() const { return params_; }

Libs/Optimize/OptimizeParameters.h

@@ -138,6 +138,8 @@ class OptimizeParameters {
   double get_shared_boundary_weight();
   void set_shared_boundary_weight(double value);
 
+  Parameters get_parameters() const;
+
  private:
   std::string get_output_prefix();
 

Python/shapeworks/shapeworks/network_analysis.py

@@ -297,7 +297,12 @@ def run(self):
         snpm = spm1d.stats.nonparam.ttest2(
             np.transpose(all_data[:, :, np.where(grouprs == 0)[0], :].reshape(num_pts, len(np.where(grouprs == 0)[0]))),
             np.transpose(all_data[:, :, np.where(grouprs == 1)[0], :].reshape(num_pts, len(np.where(grouprs == 1)[0]))))
-        snpmi = snpm.inference(0.05, two_tailed=True, iterations=n_iter, force_iterations=True)  # assume one-sided
+
+        alpha = 0.05
+        # Specified alpha must be greater than or equal to (1/nPermTotal)=0.33333
+        if (alpha < 1/n_iter):
+            alpha = 1/n_iter
+        snpmi = snpm.inference(alpha, two_tailed=True, iterations=n_iter, force_iterations=True)  # assume one-sided
 
         Z = snpmi.z  # flattened test statistic  (i.e., t value) over only non-zero-variance nodes
         tradzstar = snpmi.zstar  # critical test statistic  (i.e., critical t value)

Studio/Analysis/AnalysisTool.cpp

@@ -345,6 +345,8 @@ void AnalysisTool::set_session(QSharedPointer<Session> session) {
   ui_->group1_button->setChecked(false);
   ui_->group2_button->setChecked(false);
   update_difference_particles();
+  ui_->pcaSlider->setValue(0);
+  ui_->group_slider->setValue(10);
 
   ui_->show_predicted_scalar->setChecked(false);
   ui_->show_difference_to_predicted_scalar->setChecked(false);
@@ -380,6 +382,7 @@ void AnalysisTool::handle_analysis_options() {
     ui_->pcaAnimateCheckBox->setEnabled(false);
     ui_->pcaModeSpinBox->setEnabled(false);
     pca_animate_timer_.stop();
+    group_animate_timer_.stop();
     ui_->pcaSlider->setEnabled(false);
     if (ui_->singleSamplesRadio->isChecked()) {
       // one sample mode
@@ -406,6 +409,7 @@ void AnalysisTool::handle_analysis_options() {
     ui_->pcaSlider->setEnabled(true);
     ui_->pcaAnimateCheckBox->setEnabled(true);
     ui_->pcaModeSpinBox->setEnabled(true);
+    group_animate_timer_.stop();
     auto domain_names = session_->get_project()->get_domain_names();
     bool multiple_domains = domain_names.size() > 1;
     if (multiple_domains) {
@@ -422,6 +426,7 @@ void AnalysisTool::handle_analysis_options() {
     ui_->pcaAnimateCheckBox->setEnabled(false);
     ui_->pcaModeSpinBox->setEnabled(false);
     pca_animate_timer_.stop();
+    group_animate_timer_.stop();
   }
 
   update_difference_particles();
@@ -519,6 +524,11 @@ void AnalysisTool::group_p_values_clicked() {
 
 //---------------------------------------------------------------------------
 void AnalysisTool::network_analysis_clicked() {
+  if (ui_->network_feature->currentText().isEmpty()) {
+    QMessageBox::warning(this, "Network Analysis", "Project must have a scalar feature for network analysis.");
+    // SW_WARN("Project must have a scalar features for network analysis");
+    return;
+  }
   network_analysis_job_ =
       QSharedPointer<NetworkAnalysisJob>::create(session_->get_project(), ui_->group_combo->currentText().toStdString(),
                                                  ui_->network_feature->currentText().toStdString());
@@ -602,6 +612,7 @@ bool AnalysisTool::compute_stats() {
     if (particles.size() == 0) {
       continue;  // skip any that don't have particles
     }
+
     if (groups_active()) {
       auto group = shape->get_subject()->get_group_value(group_set);
       if (group == left_group) {
@@ -718,7 +729,6 @@ bool AnalysisTool::compute_stats() {
 //-----------------------------------------------------------------------------
 Particles AnalysisTool::get_mean_shape_points() {
   if (!compute_stats()) {
-    std::cerr << "Non buenas, returning empty particles\n";
     return Particles();
   }
 
@@ -926,6 +936,7 @@ void AnalysisTool::store_settings() {
 //---------------------------------------------------------------------------
 void AnalysisTool::shutdown() {
   pca_animate_timer_.stop();
+  group_animate_timer_.stop();
 
   for (const auto& worker : workers_) {
     if (worker) {
@@ -1453,6 +1464,8 @@ void AnalysisTool::update_group_boxes() {
 //---------------------------------------------------------------------------
 void AnalysisTool::update_group_values() {
   block_group_change_ = true;
+  stats_ready_ = false;
+
   auto values = session_->get_project()->get_group_values(ui_->group_combo->currentText().toStdString());
 
   if (values != current_group_values_) {
@@ -2107,7 +2120,9 @@ void AnalysisTool::reconstruction_method_changed() {
 void AnalysisTool::set_active(bool active) {
   if (!active) {
     ui_->pcaAnimateCheckBox->setChecked(false);
+    ui_->group_animate_checkbox->setChecked(false);
     pca_animate_timer_.stop();
+    group_animate_timer_.stop();
   } else {
     auto features = session_->get_project()->get_feature_names();
     ui_->network_feature->clear();
@@ -2137,6 +2152,12 @@ Particles AnalysisTool::convert_from_combined(const Eigen::VectorXd& points) {
   int idx = 0;
   for (int d = 0; d < worlds.size(); d++) {
     Eigen::VectorXd new_world(worlds[d].size());
+
+    if (idx + new_world.size() > points.size()) {
+      SW_WARN("Inconsistent number of values in particle vector");
+      return {};
+    }
+
     for (int i = 0; i < worlds[d].size(); i++) {
       new_world[i] = points[idx++];
     }

Studio/Groom/GroomTool.cpp

@@ -591,6 +591,18 @@ void GroomTool::on_run_groom_button_clicked() {
   timer_.start();
 
   SW_LOG("Please wait: running groom step...");
+
+  // log parameters
+  SW_LOG_ONLY("Grooming parameters:");
+
+  for (const auto& domain_name : session_->get_project()->get_domain_names()) {
+    auto params = GroomParameters(session_->get_project(), domain_name);
+    SW_LOG_ONLY("Domain: " + domain_name);
+    for (const auto& pair : params.get_parameters().get_map()) {
+      SW_LOG_ONLY(pair.first + ": " + pair.second);
+    }
+  }
+
   Q_EMIT progress(0);
 
   groom_ = QSharedPointer<Groom>(new Groom(session_->get_project()));

Studio/Interface/ShapeWorksStudioApp.cpp

@@ -442,6 +442,9 @@ void ShapeWorksStudioApp::import_files(QStringList file_names) {
   } catch (std::runtime_error& e) {
     handle_error(e.what());
   }
+
+  visualizer_->reset_camera();
+
   handle_message("Files loaded");
   handle_progress(100);
 }

Studio/Interface/UpdateChecker.cpp

@@ -50,7 +50,7 @@ void UpdateChecker::run_auto_update_check() {
 
 void UpdateChecker::run_update_check() {
   // check if a new version is available
-  QNetworkRequest request(QUrl("http://www.sci.utah.edu/~shapeworks/version.json"));
+  QNetworkRequest request(QUrl("https://www.sci.utah.edu/~shapeworks/version.json"));
   network_.get(request);
 }