Add TIME_SCOPE profiling instrumentation across major C++ pipelines

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/Application/Job/Job.cpp

@@ -1,4 +1,5 @@
 #include <Job/Job.h>
+#include <Profiling.h>
 namespace shapeworks {
 
 //---------------------------------------------------------------------------
@@ -19,6 +20,12 @@ QString Job::get_abort_message() {
   return name() + " aborted.  Duration: " + duration + " seconds";
 }
 
+//---------------------------------------------------------------------------
+void Job::execute() {
+  TIME_SCOPE(name());
+  run();
+}
+
 //---------------------------------------------------------------------------
 void Job::start_timer() { this->timer_.start(); }
 

Libs/Application/Job/Job.h

@@ -15,6 +15,9 @@ class Job : public QObject {
   //! run the job
   virtual void run() = 0;
 
+  //! execute the job with profiling instrumentation
+  void execute();
+
   //! get the name of the job
   virtual QString name() = 0;
 

Libs/Application/Job/PythonWorker.cpp

@@ -104,7 +104,7 @@ void PythonWorker::start_job(QSharedPointer<Job> job) {
       }
       Q_EMIT job->progress(0);
       current_job_ = job;
-      current_job_->run();
+      current_job_->execute();
       current_job_->set_complete(true);
       if (!job->get_quiet_mode()) {
         SW_LOG(current_job_->get_completion_message().toStdString());

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();