Add TIME_SCOPE profiling instrumentation across groom, mesh, image, analysis, and CLI pipelines

Author: akenmorris

Applications/shapeworks/Command.cpp

@@ -1,4 +1,5 @@
 #include "Command.h"
+#include <Profiling.h>
 #include <sstream>
 
 namespace shapeworks {
@@ -22,6 +23,7 @@ std::vector<std::string> Command::parse_args(const std::vector<std::string> &arg
 ///////////////////////////////////////////////////////////////////////////////
 int Command::run(SharedCommandData &sharedData)
 {
+  TIME_SCOPE(QString::fromStdString(name()));
   const optparse::Values &options = parser.get_parsed_options();
 
   return this->execute(options, sharedData) ? EXIT_SUCCESS : EXIT_FAILURE;

Libs/Alignment/Procrustes3D.cpp

@@ -1,5 +1,6 @@
 #include "Procrustes3D.h"
 
+#include <Profiling.h>
 #include <vnl/algo/vnl_svd.h>
 
 #include <iostream>
@@ -44,6 +45,7 @@ void Procrustes3D::RemoveTranslation(SimilarityTransformListType& transforms, Sh
 
 //---------------------------------------------------------------------------
 void Procrustes3D::AlignShapes(SimilarityTransformListType& transforms, ShapeListType& shapes) {
+  TIME_SCOPE("Procrustes3D::AlignShapes");
   const RealType SOS_EPSILON = 1.0e-8;
 
   PointType center;

Libs/Analyze/Analyze.cpp

@@ -2,6 +2,7 @@
 
 #include <Groom/GroomParameters.h>
 #include <Logging.h>
+#include <Profiling.h>
 #include <MeshWarper.h>
 #include <Particles/ParticleNormalEvaluation.h>
 #include <StringUtils.h>
@@ -167,6 +168,7 @@ Analyze::Analyze(ProjectHandle project) : project_(project), mesh_manager_(new M
 
 //---------------------------------------------------------------------------
 void Analyze::run_offline_analysis(std::string outfile, float range, float steps) {
+  TIME_SCOPE("Analyze::run_offline_analysis");
   SW_LOG("ShapeWorks Offline Analysis");
   if (!project_->get_particles_present()) {
     throw std::runtime_error("Project has not been optimized, please run optimize first");
@@ -478,6 +480,7 @@ bool Analyze::update_shapes() {
 
 //---------------------------------------------------------------------------
 bool Analyze::compute_stats() {
+  TIME_SCOPE("Analyze::compute_stats");
   if (stats_ready_) {
     return true;
   }

Libs/Analyze/Reconstruction.cpp

@@ -23,6 +23,7 @@
 #include "itkNrrdImageIOFactory.h"
 #include "itkMetaImageIOFactory.h"
 #include "Reconstruction.h"
+#include <Profiling.h>
 
 #include <vtkLoopSubdivisionFilter.h>
 #include <vtkButterflySubdivisionFilter.h>
@@ -67,6 +68,7 @@ vtkSmartPointer<vtkPolyData> Reconstruction<TTransformType,TInterpolatorType, TC
         std::vector< PointArrayType > local_pts,
         std::vector< PointArrayType > global_pts,
         std::vector<std::string> distance_transform) {
+    TIME_SCOPE("Reconstruction::getDenseMean");
     if (!this->denseDone_ || !local_pts.empty() ||
             !distance_transform.empty() || !global_pts.empty()) {
         this->denseDone_ = false;
@@ -439,6 +441,7 @@ void Reconstruction<TTransformType,TInterpolatorType, TCoordRep, PixelType, Imag
         std::vector< PointArrayType > local_pts,
         std::vector< PointArrayType > global_pts,
         std::vector<std::string> distance_transform) {
+    TIME_SCOPE("Reconstruction::computeDenseMean");
     try {
         //turn the sets of global points to one sparse global mean.
         float init[] = { 0.f,0.f,0.f };

Libs/Groom/Groom.cpp

@@ -2,6 +2,7 @@
 #include <GroomParameters.h>
 #include <Image/Image.h>
 #include <Logging.h>
+#include <Profiling.h>
 #include <Mesh/Mesh.h>
 #include <Mesh/MeshUtils.h>
 #include <Optimize/Constraints/Constraints.h>
@@ -32,6 +33,7 @@ Groom::Groom(ProjectHandle project) : project_{project} {}
 
 //---------------------------------------------------------------------------
 bool Groom::run() {
+  TIME_SCOPE("Groom::run");
   ShapeWorksUtils::setup_threads();
 
   used_names_.clear();
@@ -124,6 +126,7 @@ bool Groom::run() {
 
 //---------------------------------------------------------------------------
 bool Groom::image_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
+  TIME_SCOPE("Groom::image_pipeline");
   // grab parameters
   auto params = GroomParameters(project_, project_->get_domain_names()[domain]);
 
@@ -321,6 +324,7 @@ bool Groom::run_image_pipeline(Image& image, GroomParameters params) {
 
 //---------------------------------------------------------------------------
 bool Groom::mesh_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
+  TIME_SCOPE("Groom::mesh_pipeline");
   // grab parameters
   auto params = GroomParameters(project_, project_->get_domain_names()[domain]);
 
@@ -557,6 +561,7 @@ bool Groom::get_aborted() { return abort_; }
 
 //---------------------------------------------------------------------------
 bool Groom::run_alignment() {
+  TIME_SCOPE("Groom::run_alignment");
   size_t num_domains = project_->get_number_of_domains_per_subject();
   SW_DEBUG("Running alignment, number of domains = {}", num_domains);
   auto subjects = project_->get_subjects();

Libs/Image/Image.cpp

@@ -1,6 +1,7 @@
 #include "Image.h"
 
 #include <Logging.h>
+#include <Profiling.h>
 #include <itkAntiAliasBinaryImageFilter.h>
 #include <itkBinaryFillholeImageFilter.h>
 #include <itkBinaryThresholdImageFilter.h>
@@ -104,6 +105,7 @@ Image& Image::operator=(Image&& img) {
 }
 
 Image::ImageType::Pointer Image::read(const std::string& pathname) {
+  TIME_SCOPE("Image::read");
   ImageUtils::register_itk_factories();
 
   if (pathname.empty()) {
@@ -327,6 +329,7 @@ Image::ImageType::Pointer Image::readDICOMImage(const std::string& pathname) {
 }
 
 Image& Image::write(const std::string& filename, bool compressed) {
+  TIME_SCOPE("Image::write");
   if (!this->itk_image_) {
     throw std::invalid_argument("Image invalid");
   }
@@ -351,6 +354,7 @@ Image& Image::write(const std::string& filename, bool compressed) {
 }
 
 Image& Image::antialias(unsigned iterations, double maxRMSErr, int layers) {
+  TIME_SCOPE("Image::antialias");
   if (layers < 0) {
     throw std::invalid_argument("layers must be >= 0");
   }
@@ -405,6 +409,7 @@ Image& Image::resample(const TransformPtr transform, const Point3 origin, Dims d
 }
 
 Image& Image::resample(const Vector3& spacing, Image::InterpolationType interp) {
+  TIME_SCOPE("Image::resample");
   // compute logical dimensions that keep all image data for this spacing
   Dims inputDims(this->dims());
   Vector3 inputSpacing(this->spacing());
@@ -518,6 +523,7 @@ Image& Image::pad(IndexRegion& region, PixelType value) {
 }
 
 Image& Image::pad(Dims lowerExtendRegion, Dims upperExtendRegion, PixelType value) {
+  TIME_SCOPE("Image::pad");
   using FilterType = itk::ConstantPadImageFilter<ImageType, ImageType>;
   FilterType::Pointer filter = FilterType::New();
 
@@ -632,6 +638,7 @@ Image& Image::binarize(PixelType minVal, PixelType maxVal, PixelType innerVal, P
 }
 
 Image& Image::computeDT(PixelType isoValue) {
+  TIME_SCOPE("Image::computeDT");
   using FilterType = itk::ReinitializeLevelSetImageFilter<ImageType>;
   FilterType::Pointer filter = FilterType::New();
 
@@ -732,6 +739,7 @@ Image& Image::applyIntensityFilter(double minVal, double maxVal) {
 }
 
 Image& Image::gaussianBlur(double sigma) {
+  TIME_SCOPE("Image::gaussianBlur");
   using BlurType = itk::DiscreteGaussianImageFilter<ImageType, ImageType>;
   BlurType::Pointer blur = BlurType::New();
 
@@ -744,6 +752,7 @@ Image& Image::gaussianBlur(double sigma) {
 }
 
 Image& Image::crop(PhysicalRegion region, const int padding) {
+  TIME_SCOPE("Image::crop");
   region.shrink(physicalBoundingBox());  // clip region to fit inside image
   if (!region.valid()) {
     throw std::invalid_argument("Invalid region specified (it may lie outside physical bounds of image).");
@@ -862,6 +871,7 @@ Image& Image::setCoordsys(ImageType::DirectionType coordsys) {
 }
 
 Image& Image::isolate() {
+  TIME_SCOPE("Image::isolate");
   typedef itk::Image<unsigned char, 3> IsolateType;
   typedef itk::CastImageFilter<ImageType, IsolateType> ToIntType;
   ToIntType::Pointer filter = ToIntType::New();

Libs/Mesh/Mesh.cpp

@@ -63,6 +63,7 @@
 #include "Image.h"
 #include "Libs/Optimize/Domain/Surface.h"
 #include "Logging.h"
+#include "Profiling.h"
 #include "MeshComputeThickness.h"
 #include "MeshUtils.h"
 #include "PreviewMeshQC/FEAreaCoverage.h"
@@ -273,6 +274,7 @@ Mesh& Mesh::coverage(const Mesh& otherMesh, bool allowBackIntersections, double
 }
 
 Mesh& Mesh::smooth(int iterations, double relaxation) {
+  TIME_SCOPE("Mesh::smooth");
   auto smoother = vtkSmartPointer<vtkSmoothPolyDataFilter>::New();
 
   smoother->SetInputData(this->poly_data_);
@@ -293,6 +295,7 @@ Mesh& Mesh::smooth(int iterations, double relaxation) {
 }
 
 Mesh& Mesh::smoothSinc(int iterations, double passband) {
+  TIME_SCOPE("Mesh::smoothSinc");
   auto smoother = vtkSmartPointer<vtkWindowedSincPolyDataFilter>::New();
   smoother->SetInputData(this->poly_data_);
   // minimum of 2.  See docs of vtkWindowedSincPolyDataFilter for explanation
@@ -316,6 +319,7 @@ Mesh& Mesh::smoothSinc(int iterations, double passband) {
 }
 
 Mesh& Mesh::remesh(int numVertices, double adaptivity) {
+  TIME_SCOPE("Mesh::remesh");
   // ACVD is very noisy to std::cout, even with console output set to zero
   // setting the failbit on std::cout will silence this until it's cleared below
   //  std::cout.setstate(std::ios_base::failbit);
@@ -385,6 +389,7 @@ Mesh& Mesh::reflect(const Axis& axis, const Vector3& origin) {
 }
 
 MeshTransform Mesh::createTransform(const Mesh& target, Mesh::AlignmentType align, unsigned iterations) {
+  TIME_SCOPE("Mesh::createTransform");
   return createRegistrationTransform(target, align, iterations);
 }
 
@@ -413,6 +418,7 @@ Mesh& Mesh::rotate(const double angle, const Axis axis) {
 }
 
 Mesh& Mesh::fillHoles(double hole_size) {
+  TIME_SCOPE("Mesh::fillHoles");
   auto filter = vtkSmartPointer<vtkFillHolesFilter>::New();
   filter->SetInputData(this->poly_data_);
   filter->SetHoleSize(hole_size);
@@ -430,6 +436,7 @@ Mesh& Mesh::fillHoles(double hole_size) {
 }
 
 Mesh& Mesh::clean() {
+  TIME_SCOPE("Mesh::clean");
   auto clean = vtkSmartPointer<vtkCleanPolyData>::New();
   clean->ConvertPolysToLinesOff();
   clean->ConvertLinesToPointsOff();
@@ -669,6 +676,7 @@ bool Mesh::detectTriangular() {
 }
 
 std::vector<Field> Mesh::distance(const Mesh& target, const DistanceMethod method) const {
+  TIME_SCOPE("Mesh::distance");
   if (target.numPoints() == 0 || numPoints() == 0) {
     throw std::invalid_argument("meshes must have points");
   }
@@ -832,6 +840,7 @@ bool Mesh::isPointInside(const Point3 point) const {
 }
 
 double Mesh::geodesicDistance(int source, int target) const {
+  TIME_SCOPE("Mesh::geodesicDistance");
   if (source < 0 || target < 0 || numPoints() < source || numPoints() < target) {
     throw std::invalid_argument("requested point ids outside range of points available in mesh");
   }
@@ -841,6 +850,7 @@ double Mesh::geodesicDistance(int source, int target) const {
 }
 
 Field Mesh::geodesicDistance(const Point3 landmark) const {
+  TIME_SCOPE("Mesh::geodesicDistance");
   auto distance = vtkSmartPointer<vtkDoubleArray>::New();
   distance->SetNumberOfComponents(1);
   distance->SetNumberOfTuples(numPoints());
@@ -856,6 +866,7 @@ Field Mesh::geodesicDistance(const Point3 landmark) const {
 }
 
 Field Mesh::geodesicDistance(const std::vector<Point3> curve) const {
+  TIME_SCOPE("Mesh::geodesicDistance");
   auto minDistance = vtkSmartPointer<vtkDoubleArray>::New();
   minDistance->SetNumberOfComponents(1);
   minDistance->SetNumberOfTuples(numPoints());
@@ -875,6 +886,7 @@ Field Mesh::geodesicDistance(const std::vector<Point3> curve) const {
 }
 
 Field Mesh::curvature(const CurvatureType type) const {
+  TIME_SCOPE("Mesh::curvature");
   Eigen::MatrixXd V = points();
   Eigen::MatrixXi F = faces();
 

Libs/Particles/ShapeEvaluation.cpp

@@ -1,6 +1,7 @@
 #include "ShapeEvaluation.h"
 
 #include <Logging.h>
+#include <Profiling.h>
 #include <tbb/parallel_for.h>
 
 #include <Eigen/Core>
@@ -34,6 +35,7 @@ double ShapeEvaluation::compute_compactness(const ParticleSystemEvaluation& part
 //---------------------------------------------------------------------------
 Eigen::VectorXd ShapeEvaluation::compute_full_compactness(const ParticleSystemEvaluation& particle_system,
                                                           std::function<void(float)> progress_callback) {
+  TIME_SCOPE("ShapeEvaluation::compute_full_compactness");
   const int n = particle_system.num_samples();
   const int d = particle_system.num_dims();
   const int num_modes = n - 1;  // the number of modes is one less than the number of samples
@@ -64,6 +66,7 @@ Eigen::VectorXd ShapeEvaluation::compute_full_compactness(const ParticleSystemEv
 //---------------------------------------------------------------------------
 double ShapeEvaluation::compute_generalization(const ParticleSystemEvaluation& particle_system, const int num_modes,
                                                const std::string& save_to, bool surface_distance_mode) {
+  TIME_SCOPE("ShapeEvaluation::compute_generalization");
   const long n = particle_system.num_samples();
   const long d = particle_system.num_dims();
   const Eigen::MatrixXd& p = particle_system.get_matrix();
@@ -155,6 +158,7 @@ Eigen::VectorXd ShapeEvaluation::compute_full_generalization(const ParticleSyste
                                                              std::function<void(float)> progress_callback,
                                                              std::function<bool()> check_abort,
                                                              bool surface_distance_mode) {
+  TIME_SCOPE("ShapeEvaluation::compute_full_generalization");
   const long n = particle_system.num_samples();  // number of samples
   const long d = particle_system.num_dims();     // number of dimensions (e.g. number of particles * 3)
   const Eigen::MatrixXd& p = particle_system.get_matrix();
@@ -224,6 +228,7 @@ Eigen::VectorXd ShapeEvaluation::compute_full_generalization(const ParticleSyste
 //---------------------------------------------------------------------------
 double ShapeEvaluation::compute_specificity(const ParticleSystemEvaluation& particle_system, const int num_modes,
                                             const std::string& save_to, bool surface_distance_mode) {
+  TIME_SCOPE("ShapeEvaluation::compute_specificity");
   const long n = particle_system.num_samples();
   const long d = particle_system.num_dims();
   int num_values = particle_system.get_num_values_per_particle();
@@ -327,6 +332,7 @@ Eigen::VectorXd ShapeEvaluation::compute_full_specificity(const ParticleSystemEv
                                                           std::function<void(float)> progress_callback,
                                                           std::function<bool()> check_abort,
                                                           bool surface_distance_mode) {
+  TIME_SCOPE("ShapeEvaluation::compute_full_specificity");
   const long n = particle_system.num_samples();
   const long d = particle_system.num_dims();
   const int num_values = particle_system.get_num_values_per_particle();

Studio/Data/Session.cpp

@@ -29,6 +29,7 @@
 #include <JsonProjectReader.h>
 #include <JsonProjectWriter.h>
 #include <Logging.h>
+#include <Profiling.h>
 #include <MeshManager.h>
 #include <Project/Project.h>
 #include <Shape.h>
@@ -107,6 +108,7 @@ void Session::set_parent(QWidget* parent) { parent_ = parent; }
 
 //---------------------------------------------------------------------------
 bool Session::save_project(QString filename) {
+  TIME_SCOPE("Session::save_project");
   QProgressDialog progress("Saving Project...", "Abort", 0, 100, parent_);
   progress.setWindowModality(Qt::WindowModal);
   progress.setMinimumDuration(2000);
@@ -201,6 +203,7 @@ bool Session::save_project(QString filename) {
 
 //---------------------------------------------------------------------------
 bool Session::load_project(QString filename) {
+  TIME_SCOPE("Session::load_project");
   modified_ = false;
 
 

Studio/Data/ShapeWorksWorker.cpp

@@ -2,6 +2,7 @@
 #include <Data/ShapeWorksWorker.h>
 #include <Groom/Groom.h>
 #include <Logging.h>
+#include <Profiling.h>
 #include <Optimize/Optimize.h>
 #include <Optimize/OptimizeParameters.h>
 #include <Shape.h>
@@ -31,6 +32,7 @@ void ShapeworksWorker::process() {
   switch (this->type_) {
     case ShapeworksWorker::GroomType:
       try {
+        TIME_SCOPE("studio_groom");
         this->groom_->run();
       } catch (itk::ExceptionObject& ex) {
         SW_ERROR("{}", std::string("ITK Exception: ") + ex.GetDescription());
@@ -56,6 +58,7 @@ void ShapeworksWorker::process() {
       break;
     case ShapeworksWorker::OptimizeType:
       try {
+        TIME_SCOPE("studio_optimize");
         SW_LOG("Loading data...");
         this->optimize_parameters_->set_up_optimize(this->optimize_.data());
         SW_LOG("Optimizing correspondence...");
@@ -90,6 +93,7 @@ void ShapeworksWorker::process() {
       break;
     case ShapeworksWorker::ReconstructType:
       try {
+        TIME_SCOPE("studio_reconstruct");
         SW_LOG("Warping to mean space...");
         for (int i = 0; i < this->session_->get_domains_per_shape(); i++) {
           auto shapes = this->session_->get_shapes();