diff --git a/Applications/shapeworks/Command.cpp b/Applications/shapeworks/Command.cpp
index 69c97264bd..85678f6ab7 100644
--- a/Applications/shapeworks/Command.cpp
+++ b/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;
diff --git a/Libs/Alignment/Procrustes3D.cpp b/Libs/Alignment/Procrustes3D.cpp
index dd6290c609..a23da1d060 100755
--- a/Libs/Alignment/Procrustes3D.cpp
+++ b/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;
diff --git a/Libs/Analyze/Analyze.cpp b/Libs/Analyze/Analyze.cpp
index 03f51f7126..ec32ec5ab3 100644
--- a/Libs/Analyze/Analyze.cpp
+++ b/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;
   }
diff --git a/Libs/Analyze/Reconstruction.cpp b/Libs/Analyze/Reconstruction.cpp
index 1b395617f0..422b6eeeaa 100644
--- a/Libs/Analyze/Reconstruction.cpp
+++ b/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 };
diff --git a/Libs/Application/Job/Job.cpp b/Libs/Application/Job/Job.cpp
index 469e3a7e97..2ef8e22945 100644
--- a/Libs/Application/Job/Job.cpp
+++ b/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(); }
 
diff --git a/Libs/Application/Job/Job.h b/Libs/Application/Job/Job.h
index 24bbea6ddb..c94a0a993c 100644
--- a/Libs/Application/Job/Job.h
+++ b/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;
 
diff --git a/Libs/Application/Job/PythonWorker.cpp b/Libs/Application/Job/PythonWorker.cpp
index b062a6bd3b..e8f87f18a4 100644
--- a/Libs/Application/Job/PythonWorker.cpp
+++ b/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());
diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp
index 8ce36911d5..c50e444ad5 100644
--- a/Libs/Groom/Groom.cpp
+++ b/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();
diff --git a/Libs/Image/Image.cpp b/Libs/Image/Image.cpp
index 26937abdf5..88e4e86854 100644
--- a/Libs/Image/Image.cpp
+++ b/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();
diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp
index 6023bbab74..e7c9f64203 100644
--- a/Libs/Mesh/Mesh.cpp
+++ b/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();
 
diff --git a/Libs/Particles/ShapeEvaluation.cpp b/Libs/Particles/ShapeEvaluation.cpp
index 34d96a5b99..a59572e57c 100644
--- a/Libs/Particles/ShapeEvaluation.cpp
+++ b/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();
diff --git a/Studio/Data/Session.cpp b/Studio/Data/Session.cpp
index 56726fbc10..c57295d69b 100644
--- a/Studio/Data/Session.cpp
+++ b/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;
 
 
diff --git a/Studio/Data/ShapeWorksWorker.cpp b/Studio/Data/ShapeWorksWorker.cpp
index 13a1aca8c0..564c702cd2 100644
--- a/Studio/Data/ShapeWorksWorker.cpp
+++ b/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();
diff --git a/Studio/Data/Worker.cpp b/Studio/Data/Worker.cpp
index 47672ad944..46856e1def 100644
--- a/Studio/Data/Worker.cpp
+++ b/Studio/Data/Worker.cpp
@@ -42,7 +42,7 @@ void Worker::stop() {
 void Worker::process() {
   try {
     job_->start_timer();
-    job_->run();
+    job_->execute();
     job_->set_complete(true);
   } catch (std::exception& e) {
     SW_ERROR("{}", e.what());
diff --git a/Studio/Interface/ShapeWorksStudioApp.cpp b/Studio/Interface/ShapeWorksStudioApp.cpp
index 88046a9dbc..6923ecd7b3 100644
--- a/Studio/Interface/ShapeWorksStudioApp.cpp
+++ b/Studio/Interface/ShapeWorksStudioApp.cpp
@@ -20,6 +20,7 @@
 // shapeworks
 #include <Applications/Configuration.h>
 #include <Logging.h>
+#include <Profiling.h>
 #include <Mesh/Mesh.h>
 #include <SurfaceReconstructor.h>
 #include <Utils/StringUtils.h>
@@ -1526,6 +1527,7 @@ void ShapeWorksStudioApp::on_view_mode_combobox_currentIndexChanged(QString disp
 
 //---------------------------------------------------------------------------
 void ShapeWorksStudioApp::open_project(QString filename) {
+  TIME_SCOPE("open_project");
   preferences_.set_last_directory(QFileInfo(filename).absolutePath());
   new_session();
   SW_LOG("Loading Project: " + filename.toStdString());
@@ -1932,6 +1934,7 @@ void ShapeWorksStudioApp::update_alignment_options() {
 
 //---------------------------------------------------------------------------
 void ShapeWorksStudioApp::save_project(QString filename) {
+  TIME_SCOPE("save_project");
   session_->parameters().set(ShapeWorksStudioApp::SETTING_ZOOM_C, std::to_string(ui_->zoom_slider->value()));
   session_->parameters().set("notes", data_tool_->get_notes());
   session_->parameters().set("analysis_mode", analysis_tool_->get_analysis_mode());
diff --git a/docs/about/release-notes.md b/docs/about/release-notes.md
index c872c6d505..6c2ae5d286 100644
--- a/docs/about/release-notes.md
+++ b/docs/about/release-notes.md
@@ -3,22 +3,35 @@
 
 ## ShapeWorks 6.7.0
 
+![](../img/about/release6.7.png)
+
 ### What is new?
 
+  * **DeepSSM**
+    * DeepSSM CLI: `shapeworks deepssm` command for headless/batch workflows on HPC clusters
+    * TL-DeepSSM support from CLI via `--tl_net` flag
+    * Streaming data loaders to reduce memory usage during training
+    * Config schema validation before running
+    * Improved error handling in data loaders
+
   * **ShapeWorks Back-end**
     * Fixed domain optimization speedups: shape space precomputation, Procrustes caching, and batch particle loading
     * PCA projection correspondence for fixed domain optimization (`use_pca_projection` parameter)
+    * Geodesic walk for particle splitting on meshes, keeping particles on-surface
+    * Laplacian mesh warping as alternative to biharmonic, reducing thinning artifacts on sparse regions
+    * Replaced VNL eigensystem with Eigen SVD for PCA computation, fixing eigenvector sign ambiguity
+    * Parallel mesh loading with TBB and progress reporting
     * Automatic mesh repair in grooming pipeline (replaces manual Extract Largest Component step)
     * Auto ICP alignment subset size defaults to 30 shapes for faster grooming on large datasets
     * Support for multiple shared boundaries between domains
+    * Sampling gradient scaling based on surface area normalization for multi-domain models
     * Early stopping via morphological deviation score for optimization convergence detection
-    * DeepSSM CLI: `shapeworks deepssm` command for headless/batch DeepSSM workflows
-    * Laplacian mesh warping as alternative to biharmonic, reducing thinning artifacts on sparse regions
     * Built-in performance profiling and Chrome trace output via `SW_TIME_PROFILE` and `SW_TIME_TRACE` environment variables (see [developer docs](../dev/build.md#performance-profiling))
 
   * **ShapeWorks Front-end**
     * DWD (Distance Weighted Discrimination) group analysis method alongside LDA
-    * Redesigned Optimize panel
+    * Redesigned Optimize panel with reorganized parameter groups
+    * Sampling scale UI for multi-domain particle density control
     * Export clipped meshes
     * Warp method selector (Biharmonic/Laplacian) in Analysis panel
 
@@ -29,6 +42,16 @@
     * PyTorch 2.8.0 (2.2.2 on Intel Mac)
     * Linux builds use manylinux_2_28 (GLIBC 2.28, compatible with RHEL 8+, Ubuntu 20.04+)
 
+### Fixes
+  * Fix cross-domain p-value bleed in multi-domain models where neighboring domain particles corrupted KD-tree interpolation near anatomically close regions
+  * Fix multi-domain overlap when pre-aligned data skips grooming alignment (#2514)
+  * Fix reflection being applied twice in global multi-domain ICP alignment
+  * Fix PPCA noise variance dilution and null-space leakage in MorphologicalDeviationScore (#2495)
+  * Fix fillHoles creating zero-area triangles
+  * Fix crash when geodesic distance from landmarks is enabled without landmarks
+  * Fix DeepSSM crash when training split is set to 0% (#2398)
+  * Fix crash when mesh contains no cells during grooming
+
 
 ## ShapeWorks 6.6.1 - 2025-05