From 9659d7883e3fa413d0a62562f53b340479ba04d5 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Wed, 26 Aug 2026 12:34:09 -0600 Subject: [PATCH 1/5] Keep registration based initialization from diverging on dissimilar shapes The similarity metric averages over only the samples that land inside the moving image, so pulling the two images apart improves the score until nothing but background is being compared against background. The linear stages walked into that dead zone and came back mapping every particle far outside the target, which surfaced much later as an unplaceable particle from inside the sampler. A level now stops as soon as a step throws most of the fixed image off the moving one, and the affine stage is seeded with the size difference between the two shapes, which the rigid stage before it has no freedom to take up. Particles that still land where a domain cannot be sampled are walked back onto the shape rather than aborting the run, and the error names the domain and says whether the point was outside the image or merely outside the narrow band, which call for different answers. Registration initialization moves out of Optimize into RegistrationInitializer, and the transform cache key now carries the registration's own account of its settings, so that changing the algorithm cannot leave a stale description behind in the caller. --- Libs/Image/ImageRegistration.cpp | 206 +++++++- Libs/Image/ImageRegistration.h | 11 + Libs/Optimize/Domain/ImageDomain.h | 56 +- Libs/Optimize/Domain/ParticleDomain.h | 4 + Libs/Optimize/Optimize.cpp | 420 +-------------- Libs/Optimize/Optimize.h | 33 +- Libs/Optimize/RegistrationInitializer.cpp | 596 ++++++++++++++++++++++ Libs/Optimize/RegistrationInitializer.h | 81 +++ 8 files changed, 952 insertions(+), 455 deletions(-) create mode 100644 Libs/Optimize/RegistrationInitializer.cpp create mode 100644 Libs/Optimize/RegistrationInitializer.h diff --git a/Libs/Image/ImageRegistration.cpp b/Libs/Image/ImageRegistration.cpp index cd4dcf4e46b..ef5716edf91 100644 --- a/Libs/Image/ImageRegistration.cpp +++ b/Libs/Image/ImageRegistration.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -44,26 +45,43 @@ using LinearMetricType = itk::MeanSquaresImageToImageMetricv4; using OptimizerType = itk::ConjugateGradientLineSearchOptimizerv4Template; -constexpr double kLineSearchLowerLimit = 0.0; -constexpr double kLineSearchUpperLimit = 2.0; -constexpr double kLineSearchEpsilon = 0.2; +//! Names this registration, for callers that keep transforms between runs. Change it whenever the +//! algorithm changes, to a name that has not been used before: a bare version number invites reuse +//! of one already written to somebody's cache, and a stale transform then reads back as a hit, which +//! is far harder to notice than a miss. +constexpr const char* ALGORITHM_NAME = "syn-overlap-guard-affine-scale-seed"; + +constexpr double LINE_SEARCH_LOWER_LIMIT = 0.0; +constexpr double LINE_SEARCH_UPPER_LIMIT = 2.0; +constexpr double LINE_SEARCH_EPSILON = 0.2; // A heavily shrunk level of an already small image carries almost no signal, and the optimizer // responds by taking a large wrong step that the finer levels never recover from. Never shrink a // dimension below this many voxels. -constexpr unsigned int kMinimumLevelSize = 16; +constexpr unsigned int MINIMUM_LEVEL_SIZE = 16; // The furthest a single optimizer step may move a point, as a fraction of the largest physical // dimension of the images being registered. -constexpr double kMaximumStepFraction = 0.1; +constexpr double MAXIMUM_STEP_FRACTION = 0.1; + +// How much of the overlap between the two images a step is allowed to give up before it is treated +// as having lost the shape it was registering. See OverlapGuardCommand. +constexpr double MINIMUM_OVERLAP_FRACTION = 0.5; + +// The range of size differences the affine stage will be seeded with. Shapes of the same anatomy do +// not differ by more than this, so an estimate outside it is a sign that the images are not what the +// estimate assumes -- and the seed is only ever a starting point worth having, never one worth +// insisting on. +constexpr double MINIMUM_SEED_SCALE = 0.5; +constexpr double MAXIMUM_SEED_SCALE = 2.0; //--------------------------------------------------------------------------- -/// Reduce the requested shrink factors so that no level shrinks the image below kMinimumLevelSize. +/// Reduce the requested shrink factors so that no level shrinks the image below MINIMUM_LEVEL_SIZE. std::vector clamp_shrink_factors(const std::vector& shrink_factors, const ImageType* image) { const auto size = image->GetBufferedRegion().GetSize(); const auto smallest = std::min({size[0], size[1], size[2]}); - const auto largest_useful = std::max(1u, static_cast(smallest / kMinimumLevelSize)); + const auto largest_useful = std::max(1u, static_cast(smallest / MINIMUM_LEVEL_SIZE)); std::vector clamped; clamped.reserve(shrink_factors.size()); @@ -73,6 +91,130 @@ std::vector clamp_shrink_factors(const std::vector& return clamped; } +//--------------------------------------------------------------------------- +/// How much bigger the moving shape is than the fixed one, and the point to scale about, taken from +/// the second moments of the two images. +/// +/// The rigid stage runs first and cannot scale at all, so where the two shapes differ in size the +/// best it can do is put the smaller somewhere inside the larger -- an underdetermined problem with +/// no single answer. The affine stage then starts from a pose that says nothing about the size +/// difference, and on shapes far enough apart in size it never finds it: the registration comes back +/// having mapped the reference onto part of the target and left the rest bare. The size difference +/// is the one thing here that can be measured directly rather than searched for, so measure it and +/// hand it to the affine stage as its starting point. +struct ScaleSeed { + bool usable{false}; + double scale{1.0}; + ImageType::PointType center; +}; + +ScaleSeed estimate_scale_seed(const ImageType* fixed, const ImageType* moving) { + using MomentsType = itk::ImageMomentsCalculator; + + auto measure = [](const ImageType* image) { + auto moments = MomentsType::New(); + moments->SetImage(const_cast(image)); + moments->Compute(); + return moments; + }; + + ScaleSeed seed; + auto fixed_moments = measure(fixed); + auto moving_moments = measure(moving); + + const double fixed_mass = fixed_moments->GetTotalMass(); + const double moving_mass = moving_moments->GetTotalMass(); + if (!(fixed_mass > 0.0) || !(moving_mass > 0.0)) { + return seed; + } + + // ITK returns each principal moment as the total mass times the mean squared radius about that + // axis, so dividing the mass back out leaves a length that the two images can be compared by + const auto fixed_moment = fixed_moments->GetPrincipalMoments(); + const auto moving_moment = moving_moments->GetPrincipalMoments(); + + double product = 1.0; + for (unsigned int i = 0; i < 3; i++) { + const double fixed_radius = std::sqrt(fixed_moment[i] / fixed_mass); + const double moving_radius = std::sqrt(moving_moment[i] / moving_mass); + if (!(fixed_radius > 0.0) || !(moving_radius > 0.0)) { + return seed; + } + product *= moving_radius / fixed_radius; + } + + // one scale for all three axes: pairing the axes individually would need to know which of the + // moving shape's axes answers to which of the fixed shape's, and on anything near round they + // cannot be told apart + const double scale = std::cbrt(product); + if (!std::isfinite(scale) || scale < MINIMUM_SEED_SCALE || scale > MAXIMUM_SEED_SCALE) { + return seed; + } + + const auto center_of_gravity = fixed_moments->GetCenterOfGravity(); + for (unsigned int i = 0; i < 3; i++) { + seed.center[i] = center_of_gravity[i]; + } + seed.scale = scale; + seed.usable = true; + return seed; +} + +//--------------------------------------------------------------------------- +/// Stop a level as soon as a step throws most of the fixed image off the moving one. +/// +/// The metric averages over only those samples that land inside the moving image; ones that fall +/// outside are dropped rather than penalised. Pulling the two images apart therefore *improves* the +/// score, and in the limit leaves nothing being compared but background against background, which is +/// a perfect match. That dead zone is downhill all the way from a correct alignment, so an optimizer +/// that wanders into it stays, and the registration comes back mapping every point far outside the +/// image it was supposed to land on. +/// +/// How many samples the metric still had to work with is the tell. When that collapses, stop: the +/// optimizer holds on to the best parameters it found, which are the ones from before the escape, +/// and the remaining levels carry on from there. Stopping a level that legitimately needed to give +/// up half its overlap only costs some convergence, which is much the cheaper mistake. +template +class OverlapGuardCommand : public itk::Command { + public: + using Self = OverlapGuardCommand; + using Pointer = itk::SmartPointer; + itkNewMacro(Self); + + void set_metric(const TMetric* metric) { metric_ = metric; } + + void Execute(itk::Object* caller, const itk::EventObject& event) override { + Execute(const_cast(caller), event); + } + + void Execute(const itk::Object* caller, const itk::EventObject& event) override { + auto* optimizer = const_cast(dynamic_cast(caller)); + if (!optimizer || !metric_) { + return; + } + + if (!itk::IterationEvent().CheckEvent(&event)) { + return; + } + + // the count rises as the registration moves up the pyramid, so compare against the most this + // registration has managed rather than against a fixed number + const auto valid = metric_->GetNumberOfValidPoints(); + most_valid_ = std::max(most_valid_, valid); + + if (valid < static_cast(most_valid_ * MINIMUM_OVERLAP_FRACTION)) { + optimizer->StopOptimization(); + } + } + + protected: + OverlapGuardCommand() = default; + + private: + const TMetric* metric_{nullptr}; + itk::SizeValueType most_valid_{0}; +}; + //--------------------------------------------------------------------------- /// Build an optimizer with the scale estimator wired to the given metric, so that rotation, /// translation and scaling parameters all step by a comparable physical distance. @@ -96,14 +238,19 @@ OptimizerType::Pointer make_optimizer(TMetric* metric, unsigned int iterations, auto optimizer = OptimizerType::New(); optimizer->SetNumberOfIterations(iterations); optimizer->SetScalesEstimator(scales_estimator); - optimizer->SetMaximumStepSizeInPhysicalUnits(largest_extent * kMaximumStepFraction); + optimizer->SetMaximumStepSizeInPhysicalUnits(largest_extent * MAXIMUM_STEP_FRACTION); optimizer->SetDoEstimateLearningRateOnce(false); optimizer->SetDoEstimateLearningRateAtEachIteration(true); - optimizer->SetLowerLimit(kLineSearchLowerLimit); - optimizer->SetUpperLimit(kLineSearchUpperLimit); - optimizer->SetEpsilon(kLineSearchEpsilon); + optimizer->SetLowerLimit(LINE_SEARCH_LOWER_LIMIT); + optimizer->SetUpperLimit(LINE_SEARCH_UPPER_LIMIT); + optimizer->SetEpsilon(LINE_SEARCH_EPSILON); // keep whatever scored best rather than wherever the last step happened to land optimizer->SetReturnBestParametersAndValue(true); + + auto guard = OverlapGuardCommand::New(); + guard->set_metric(metric); + optimizer->AddObserver(itk::IterationEvent(), guard); + return optimizer; } @@ -254,13 +401,23 @@ void ImageRegistration::Impl::run_affine() { auto metric = LinearMetricType::New(); using RegistrationType = itk::ImageRegistrationMethodv4; + // Start from the size difference between the two shapes rather than from the identity. The + // composite applies what it holds in reverse, so this affine acts on points still in fixed space, + // and the point to scale about is the fixed image's own centre of gravity. + auto affine = AffineTransformType::New(); + const auto seed = estimate_scale_seed(fixed, moving); + if (seed.usable) { + affine->SetCenter(seed.center); + affine->Scale(seed.scale); + } + auto registration = RegistrationType::New(); registration->SetFixedImage(fixed); registration->SetMovingImage(moving); registration->SetMetric(metric); // the rigid result carries the coarse alignment; the affine stage only has to find the residual registration->SetMovingInitialTransform(composite); - registration->SetInitialTransform(AffineTransformType::New()); + registration->SetInitialTransform(affine); registration->InPlaceOn(); registration->SetOptimizer(make_optimizer(metric.GetPointer(), linear_iterations.front(), fixed)); apply_multi_resolution_schedule(registration.GetPointer(), effective_linear_shrink_factors, linear_smoothing_sigmas); @@ -507,6 +664,31 @@ ImageRegistration::CompositeTransformType::Pointer ImageRegistration::get_transf static std::string displacement_field_path(const std::string& filename) { return filename + ".field.nrrd"; } +//--------------------------------------------------------------------------- +std::string ImageRegistration::settings_description() const { + auto list = [](const auto& values) { + std::string text; + for (const auto& value : values) { + text += (text.empty() ? "" : ",") + std::to_string(value); + } + return text; + }; + + std::string description = ALGORITHM_NAME; + description += "|transform=" + std::to_string(static_cast(impl_->transform_type)); + description += "|step=" + std::to_string(impl_->gradient_step); + description += "|update=" + std::to_string(impl_->update_field_variance); + description += "|total=" + std::to_string(impl_->total_field_variance); + description += "|iterations=" + list(impl_->iterations); + description += "|shrink=" + list(impl_->shrink_factors); + description += "|sigmas=" + list(impl_->smoothing_sigmas); + description += "|linear_iterations=" + list(impl_->linear_iterations); + description += "|linear_shrink=" + list(impl_->linear_shrink_factors); + description += "|linear_sigmas=" + list(impl_->linear_smoothing_sigmas); + description += "|radius=" + std::to_string(impl_->correlation_radius); + return description; +} + //--------------------------------------------------------------------------- void ImageRegistration::save_transform(const std::string& filename) const { if (!impl_->composite) { diff --git a/Libs/Image/ImageRegistration.h b/Libs/Image/ImageRegistration.h index 299cb4a27e0..a9ea0e4a992 100644 --- a/Libs/Image/ImageRegistration.h +++ b/Libs/Image/ImageRegistration.h @@ -91,6 +91,17 @@ class ImageRegistration { //! the composed transform from all stages, mapping fixed space to moving space CompositeTransformType::Pointer get_transform() const; + /** + * @brief Describe the algorithm and the settings a transform would be computed with. + * + * Anything that would change the transform appears here, so a caller keeping transforms between + * runs can tell whether one it saved earlier is still the transform it would get today. That + * belongs with the algorithm rather than with the caller: whoever changes how registration works + * is editing this class, and a description left behind in somebody else's file is one that gets + * forgotten, which turns a stale transform into a cache hit instead of a miss. + */ + std::string settings_description() const; + /** * @brief Write the transform computed by run() so that it can be reused. * diff --git a/Libs/Optimize/Domain/ImageDomain.h b/Libs/Optimize/Domain/ImageDomain.h index 9201c978646..a6431286c42 100644 --- a/Libs/Optimize/Domain/ImageDomain.h +++ b/Libs/Optimize/Domain/ImageDomain.h @@ -46,6 +46,7 @@ class ImageDomain : public ParticleRegionDomain { modifies the parent class LowerBound and UpperBound. */ void SetImage(ImageType* I, double narrow_band) { this->m_FixedDomain = false; + m_NarrowBand = narrow_band; // this->Modified(); openvdb::initialize(); // It is safe to initialize multiple times. @@ -142,6 +143,19 @@ class ImageDomain : public ParticleRegionDomain { } } + /** Whether the distance transform is actually kept at this location: inside the image, and near + enough to the surface to fall within the narrow band, which is all that is retained. */ + bool IsValidLocation(const PointType& p) const override { + if (!m_VDBImage) { + return true; + } + if (!this->IsInsideBuffer(p)) { + return false; + } + const auto idxCoord = this->transform()->worldToIndex(openvdb::Vec3R(p[0], p[1], p[2])); + return !m_VDBImage->tree().isValueOff(openvdb::Coord::round(idxCoord)); + } + inline double GetMaxDiameter() const override { double bestRadius = 0; double maxdim = 0; @@ -198,21 +212,51 @@ class ImageDomain : public ParticleRegionDomain { // Make sure the coordinate is part of the narrow band if (m_VDBImage->tree().isValueOff( openvdb::Coord::round(idxCoord))) { // `isValueOff` requires an integer coordinate + const std::string message = DescribeUnsampleablePoint(p); // If multiple threads crash here at the same time, the error message displayed is just "terminate called // recursively", which isn't helpful. So we std::cerr the error to make sure its printed to the console. - std::cerr << "Sampled point outside the narrow band: " << p << std::endl; - - std::ostringstream message; - message << "Attempt to sample at a point outside the narrow band: " << p - << ". Consider increasing the narrow band"; - throw std::runtime_error(message.str()); + std::cerr << message << std::endl; + throw std::runtime_error(message); } return idxCoord; } + /// Explain why a point cannot be sampled. Two quite different problems arrive here and they need + /// different answers: a particle that has drifted a little way off the surface wants a wider narrow + /// band, but one that has been clamped against the edge of the image started out somewhere else + /// entirely, and no width of band would have caught it. + std::string DescribeUnsampleablePoint(const PointType& p) const { + std::ostringstream message; + message << "Domain " << m_DomainID << " (" << m_DomainName << "): cannot sample the distance transform at " << p + << ". "; + + // a particle is clamped to the bounding box before it is sampled, so one that arrived from + // outside the image is now sitting exactly on one of its faces + const double tolerance = GetSpacing().GetVnlVector().max_value(); + bool on_edge = false; + for (unsigned int i = 0; i < DIMENSION; i++) { + on_edge = on_edge || p[i] <= GetLowerBound()[i] + tolerance || p[i] >= GetUpperBound()[i] - tolerance; + } + + if (on_edge) { + message << "It lies on the edge of the image (" << GetLowerBound() << " to " << GetUpperBound() + << "), which is where a particle ends up when it is placed outside the image altogether, nowhere near " + "the surface. Widening the narrow band will not help. The particle was put in the wrong place to " + "begin with, usually by an initialization that did not converge or by shapes that grooming left " + "unaligned."; + } else { + message << "It is inside the image but further than the narrow band (" << m_NarrowBand + << ") from the surface, so no distance transform is kept there. Consider increasing the narrow band " + "optimization parameter."; + } + + return message.str(); + } + private: openvdb::FloatGrid::Ptr m_VDBImage; + double m_NarrowBand{0.0}; typename ImageType::SizeType m_Size; typename ImageType::SpacingType m_Spacing; PointType m_Origin; diff --git a/Libs/Optimize/Domain/ParticleDomain.h b/Libs/Optimize/Domain/ParticleDomain.h index 71263486e34..83c5784ffe6 100644 --- a/Libs/Optimize/Domain/ParticleDomain.h +++ b/Libs/Optimize/Domain/ParticleDomain.h @@ -72,6 +72,10 @@ class ParticleDomain { * PowerOfTwoPointTree. */ virtual const PointType &GetUpperBound() const = 0; + /** Whether a particle placed here could be sampled by this domain. Domains that can be sampled + anywhere have nothing to check and keep the default. */ + virtual bool IsValidLocation(const PointType &p) const { return true; } + /** Get any valid point on the domain. This is used to place the first particle. */ virtual PointType GetZeroCrossingPoint() const = 0; /** Use for neighborhood radius. */ diff --git a/Libs/Optimize/Optimize.cpp b/Libs/Optimize/Optimize.cpp index 2aaf26786e4..f28da76e98f 100644 --- a/Libs/Optimize/Optimize.cpp +++ b/Libs/Optimize/Optimize.cpp @@ -2,14 +2,14 @@ #include #include #include +#include #include // #include +#include #include #include #include -#include - #include "Profiling.h" #ifdef _WIN32 @@ -33,12 +33,14 @@ #include "EarlyStoppingConfig.h" #include "Libs/Optimize/Domain/MeshDomain.h" #include "Libs/Optimize/Domain/Surface.h" -#include "MeshUtils.h" #include "Libs/Optimize/Utils/ObjectReader.h" #include "Libs/Optimize/Utils/ObjectWriter.h" #include "Libs/Optimize/Utils/ParticleGoodBadAssessment.h" #include "Logging.h" +#include "MeshUtils.h" #include "Optimize.h" + +#include "RegistrationInitializer.h" #include "OptimizeParameterFile.h" #include "OptimizeParameters.h" #include "ShapeworksUtils.h" @@ -50,29 +52,9 @@ namespace py = pybind11; namespace shapeworks { -//! Width of the retained band around the surface, as a fraction of the shape's largest dimension, -//! when it is not set explicitly. The band is what the similarity metric can actually see, so it -//! must be a physical width independent of the rasterization grid: tying it to the grid would -//! silently narrow it as the grid is raised and starve the metric of context on dissimilar shapes. -//! ~5% of the shape size was the best compromise across similar and dissimilar pairs in testing -- -//! wide enough to bridge shape differences, not so wide that it loses surface detail. -static constexpr double kRegistrationBandFraction = 0.05; - -//! Width of the image-domain band, in voxels of the (fixed) distance transform. Image inputs already -//! carry their own resolution, so there is no grid to decouple from. -static constexpr double kRegistrationBandVoxels = 4.0; - -//! A transferred particle further than this fraction of the shape's size from the surface is -//! considered mislanded, whatever units the data is in -static constexpr double kTransferFarFraction = 0.05; - -//! Warn that a registration may have failed once this fraction of a shape's particles are mislanded. -//! Keyed off a fraction rather than the single worst particle so a stray outlier is not a false alarm. -static constexpr double kTransferFarFractionThreshold = 0.25; - //! Share of the reported progress reserved for the registrations, which dominate the wall clock in //! registration based initialization but run no particle iterations of their own -static constexpr double kTransferProgressShare = 0.85; +static constexpr double TRANSFER_PROGRESS_SHARE = 0.85; #ifdef _WIN32 static std::string find_in_path(std::string file) { @@ -745,390 +727,7 @@ void Optimize::Initialize() { } //--------------------------------------------------------------------------- -Mesh Optimize::GetDomainSurface(int domain) { - auto* particle_domain = m_sampler->GetParticleSystem()->GetDomain(domain); - - if (auto* mesh_domain = dynamic_cast(particle_domain)) { - return Mesh(mesh_domain->get_surface()->get_polydata()); - } - - // an image domain keeps only a narrow band of its distance transform, so read the groomed input - // back from disk to recover the full surface - if (domain >= static_cast(m_domain_paths.size()) || m_domain_paths[domain].empty()) { - throw std::runtime_error( - "Registration based initialization needs the path of each groomed input, but none was given for domain " + - std::to_string(domain)); - } - return Image(m_domain_paths[domain]).toMesh(0.0); -} - -//--------------------------------------------------------------------------- -Image Optimize::GetRegistrationImage(int domain) { - auto* particle_domain = m_sampler->GetParticleSystem()->GetDomain(domain); - - if (particle_domain->GetDomainType() == DomainType::Image) { - if (domain >= static_cast(m_domain_paths.size()) || m_domain_paths[domain].empty()) { - throw std::runtime_error( - "Registration based initialization needs the path of each groomed input, but none was given for domain " + - std::to_string(domain)); - } - // the groomed input is already a distance transform, so take its resolution as given - Image dt(m_domain_paths[domain]); - const auto spacing = dt.spacing(); - const double largest_spacing = std::max({spacing[0], spacing[1], spacing[2]}); - return ImageRegistration::make_registration_image(dt, GetRegistrationBand(largest_spacing)); - } - - // A mesh carries no resolution of its own, so pick one from how large it actually is. Assuming a - // spacing (and a band) in millimeters would break on data stored in other units, where a fixed - // spacing yields either an unusable grid or one too coarse to register. - Mesh mesh = GetDomainSurface(domain); - auto region = mesh.boundingBox(); - const auto extent = region.size(); - const double largest_extent = std::max({extent[0], extent[1], extent[2]}); - const double spacing = largest_extent / m_registration_grid_size; - - // Band is a physical width derived from the shape's size, independent of the chosen grid, so the - // grid controls only resolution and the two do not confound each other. - const double band = m_registration_band > 0.0 ? m_registration_band : kRegistrationBandFraction * largest_extent; - - // pad far enough that the whole band around the surface stays inside the grid - region.pad(band * 2.0); - return ImageRegistration::make_registration_image(mesh.toDistanceTransform(region, Point3({spacing, spacing, spacing})), - band); -} - -//--------------------------------------------------------------------------- -double Optimize::GetRegistrationBand(double spacing) const { - if (m_registration_band > 0.0) { - return m_registration_band; - } - // default to a band a few voxels wide, which keeps the metric focused near the surface without - // making it so thin that it spans less than a voxel - return kRegistrationBandVoxels * spacing; -} - -//--------------------------------------------------------------------------- -int Optimize::ResolveRegistrationReference() { - const int num_subjects = GetNumberOfSubjects(); - if (num_subjects < 1) { - throw std::runtime_error("No shapes to initialize"); - } - - if (m_registration_reference >= 0) { - if (m_registration_reference >= num_subjects) { - throw std::runtime_error("Requested registration reference " + std::to_string(m_registration_reference) + - " is out of range, there are only " + std::to_string(num_subjects) + " shapes"); - } - return m_registration_reference; - } - - if (num_subjects == 1) { - return 0; - } - - // combine each subject's domains into one mesh so that the template is representative of the whole - // anatomy rather than of a single domain - PrintStartMessage("Choosing a registration reference..."); - std::vector meshes; - for (int s = 0; s < num_subjects; s++) { - Mesh mesh = GetDomainSurface(s * m_domains_per_shape); - for (int d = 1; d < static_cast(m_domains_per_shape); d++) { - mesh += GetDomainSurface(s * m_domains_per_shape + d); - } - meshes.push_back(mesh); - } - - const int reference = MeshUtils::findReferenceMesh(meshes); - if (reference < 0 || reference >= num_subjects) { - throw std::runtime_error("Could not choose a registration reference"); - } - PrintDoneMessage(); - return reference; -} - -//--------------------------------------------------------------------------- -void Optimize::SpreadParticlesOnReference(int reference_shape) { - auto* system = m_sampler->GetParticleSystem(); - const int first_domain = reference_shape * m_domains_per_shape; - - // only the reference carries particles for the moment, so the correspondence matrices cannot keep - // a consistent layout; they are brought back once every shape has been populated - m_sampler->SetCorrespondenceMatricesSuspended(true); - - // seed a single particle on each of the reference's domains - for (int d = 0; d < static_cast(m_domains_per_shape); d++) { - const int domain = first_domain + d; - if (system->GetNumberOfParticles(domain) == 0) { - auto* particle_domain = system->GetDomain(domain); - system->AddPosition(particle_domain->GetValidLocationNear(particle_domain->GetZeroCrossingPoint()), domain); - } - } - system->SynchronizePositions(); - - // allocate everything without optimizing, so that turning correspondence off below sticks; Execute - // forces it back on the first time it runs - m_sampler->Initialize(); - - // only one shape carries particles at this point, so there is no correspondence to establish yet - m_sampler->SetCorrespondenceOff(); - - // report progress against the reference rather than the first shape, which is still empty - m_progress_domain_offset = first_domain; - - const double epsilon = m_spacing; - - auto needs_split = [&]() { - for (int d = 0; d < static_cast(m_domains_per_shape); d++) { - if (system->GetNumberOfParticles(first_domain + d) < m_number_of_particles[d]) { - return true; - } - } - return false; - }; - - while (needs_split() && !m_aborted) { - OptimizerStop(); - - for (int d = 0; d < static_cast(m_domains_per_shape); d++) { - const int domain = first_domain + d; - if (system->GetNumberOfParticles(domain) < m_number_of_particles[d]) { - system->SplitAllParticlesInDomain(epsilon, domain); - } - } - system->SynchronizePositions(); - - m_split_number++; - if (m_verbosity_level > 0) { - std::string counts; - for (int d = 0; d < static_cast(m_domains_per_shape); d++) { - counts += " " + std::to_string(system->GetNumberOfParticles(first_domain + d)); - } - SW_LOG("Reference split {}, particle count:{}", m_split_number, counts); - } - - m_energy_a.clear(); - m_energy_b.clear(); - m_total_energy.clear(); - m_str_energy = "split" + std::to_string(m_split_number) + "pts_init"; - - m_sampler->GetOptimizer()->set_maximum_number_of_iterations(m_iterations_per_split); - m_sampler->GetOptimizer()->set_number_of_iterations(0); - m_sampler->Execute(); - } - - m_progress_domain_offset = 0; -} - -//--------------------------------------------------------------------------- -std::string Optimize::GetRegistrationCachePath(int reference_domain, int domain) const { - if (m_registration_cache_dir.empty()) { - return ""; - } - - auto describe = [&](int d) { - std::string description = d < static_cast(m_domain_paths.size()) ? m_domain_paths[d] : std::to_string(d); - // a groomed input that has been rewritten must not read back as a hit - if (ShapeWorksUtils::file_exists(description)) { - try { - description += ":" + std::to_string(boost::filesystem::file_size(description)); - description += ":" + std::to_string(boost::filesystem::last_write_time(description)); - } catch (const std::exception&) { - } - } - return description; - }; - - // everything the transform depends on, and nothing that it does not - std::string key = describe(reference_domain) + "->" + describe(domain); - key += "|transform=" + std::to_string(static_cast(m_registration_transform_type)); - key += "|step=" + std::to_string(m_registration_gradient_step); - key += "|sigma=" + std::to_string(m_registration_flow_sigma); - key += "|band=" + std::to_string(m_registration_band); - key += "|grid=" + std::to_string(m_registration_grid_size); - - const auto hash = std::hash{}(key); - std::stringstream name; - name << m_registration_cache_dir << "/registration_" << std::hex << hash << ".tfm"; - return name.str(); -} - -//--------------------------------------------------------------------------- -void Optimize::TransferParticlesFromReference(int reference_shape) { - auto* system = m_sampler->GetParticleSystem(); - const int num_subjects = GetNumberOfSubjects(); - - for (int d = 0; d < static_cast(m_domains_per_shape) && !m_aborted; d++) { - const int reference_domain = reference_shape * m_domains_per_shape + d; - - std::vector reference_points; - for (auto k = 0; k < system->GetNumberOfParticles(reference_domain); k++) { - reference_points.push_back(system->GetPosition(k, reference_domain)); - } - - // The reference image is only needed to run a registration, so it is built lazily on the first - // cache miss. Rasterizing a large mesh takes several seconds, and when every transfer is a cache - // hit (e.g. re-running with different optimization parameters) it is never needed at all. - std::optional reference_image; - - for (int s = 0; s < num_subjects && !m_aborted; s++) { - if (s == reference_shape) { - continue; - } - const int domain = s * m_domains_per_shape + d; - - UpdateProgress(fmt::format("Registering shape {} of {}", s + 1, num_subjects)); - RefreshDuringTransfer(); - - ImageRegistration registration; - registration.set_transform_type(m_registration_transform_type); - registration.set_gradient_step(m_registration_gradient_step); - registration.set_update_field_variance(m_registration_flow_sigma); - - const auto cache_path = GetRegistrationCachePath(reference_domain, domain); - // check for the file first: asking the reader for one that is not there is a normal cache - // miss, but it makes the HDF5 layer print an alarming block of text - const bool cached = - !cache_path.empty() && ShapeWorksUtils::file_exists(cache_path) && registration.load_transform(cache_path); - if (cached) { - SW_DEBUG("Reusing cached registration: {}", cache_path); - } else { - if (!reference_image) { - // only now, on a genuine miss, is the reference image worth building - UpdateProgress(m_domains_per_shape > 1 ? fmt::format("Preparing registration (domain {})", d + 1) - : std::string("Preparing registration")); - RefreshDuringTransfer(); - reference_image = GetRegistrationImage(reference_domain); - } - registration.run(*reference_image, GetRegistrationImage(domain)); - if (!cache_path.empty()) { - try { - registration.save_transform(cache_path); - } catch (const std::exception& e) { - SW_WARN("Unable to cache registration: {}", e.what()); - } - } - } - - auto transferred = registration.transform_points(reference_points); - - // AddPositionList applies the domain constraints, which pulls each point onto the surface - system->AddPositionList(transferred, domain); - - ReportTransferQuality(domain, reference_points, transferred); - - // registrations do not run particle iterations, so charge each one its share of the budget - // reserved for them; without this the bar would sit still through the longest phase - current_particle_iterations_ += m_transfer_iteration_weight; - UpdateProgress(fmt::format("Registered shape {} of {}", s + 1, num_subjects)); - // the shape now has its particles; let the GUI redraw them (and the status) between - // registrations rather than only when the whole phase ends - RefreshDuringTransfer(); - } - } -} - -//--------------------------------------------------------------------------- -void Optimize::ReportTransferQuality(int domain, const std::vector& reference_points, - const std::vector& transferred) { - if (transferred.empty()) { - return; - } - - auto* system = m_sampler->GetParticleSystem(); - - // Judge the transfer against the size of the shape rather than an absolute distance, so the same - // thresholds work whatever units and resolution the data is in. Take that size from the domain - // itself, not from the particles: particles that have collapsed together would otherwise shrink - // the scale in step with the distances being judged, and always look acceptable. - const auto* particle_domain = system->GetDomain(domain); - const double shape_scale = particle_domain->GetLowerBound().EuclideanDistanceTo(particle_domain->GetUpperBound()); - - // A particle beyond this fraction of the shape's size is clearly in the wrong place - const double far_distance = kTransferFarFraction * shape_scale; - - double total_snap = 0.0; - double worst_snap = 0.0; - int far_count = 0; - int unmoved = 0; - - for (size_t i = 0; i < transferred.size(); i++) { - // how far the point had to travel to reach the surface is how far off the surface it landed - const double snap = transferred[i].EuclideanDistanceTo(system->GetPosition(i, domain)); - total_snap += snap; - worst_snap = std::max(worst_snap, snap); - if (shape_scale > 0.0 && snap > far_distance) { - far_count++; - } - - // a point outside the displacement field is returned unchanged rather than reported as an error, - // so an unmoved point means the field did not cover it - if (transferred[i] == reference_points[i]) { - unmoved++; - } - } - - const double mean_snap = total_snap / transferred.size(); - const double far_fraction = static_cast(far_count) / transferred.size(); - // a registration that failed leaves many particles far from the surface, not just one outlier - const bool suspect = far_fraction > kTransferFarFractionThreshold || unmoved > 0; - - const std::string name = domain < static_cast(m_filenames.size()) ? m_filenames[domain] : std::to_string(domain); - - if (m_verbosity_level > 0 || suspect) { - SW_LOG("{}: transferred particles landed {:.3g} from the surface on average (worst {:.3g})", name, mean_snap, - worst_snap); - } - - if (unmoved > 0) { - SW_WARN("{}: {} of {} transferred particles fell outside the registration field and did not move", name, unmoved, - transferred.size()); - } - - if (far_fraction > kTransferFarFractionThreshold) { - SW_WARN("{}: {:.0f}% of transferred particles landed far from the surface, the registration may have failed", name, - 100.0 * far_fraction); - } -} - -//--------------------------------------------------------------------------- -void Optimize::InitializeFromRegistration() { - m_registration_reference_chosen = ResolveRegistrationReference(); - - const std::string name = m_registration_reference_chosen * m_domains_per_shape < m_filenames.size() - ? m_filenames[m_registration_reference_chosen * m_domains_per_shape] - : std::to_string(m_registration_reference_chosen); - SW_LOG("Spreading particles on reference shape {} ({})", m_registration_reference_chosen, name); - - SpreadParticlesOnReference(m_registration_reference_chosen); - - // each shape is populated with a full set below, so the matrices can track them again - m_sampler->SetCorrespondenceMatricesSuspended(false); - - if (!m_aborted) { - TransferParticlesFromReference(m_registration_reference_chosen); - } - - // Only now does every shape hold its particles. The matrices size themselves from the first - // shape's domains, so they cannot be brought up to date any earlier: while particles were being - // spread, only the reference (which is usually not the first shape) had any. - auto* system = m_sampler->GetParticleSystem(); - system->ResyncObservers(); - system->SynchronizePositions(); - - // every shape now carries a full set of corresponding particles - m_sampler->SetCorrespondenceOn(); - - this->WritePointFiles(); - this->WritePointFilesWithFeatures(); - this->WriteTransformFile(); - this->WriteTransformFiles(); - this->WriteCuttingPlanePoints(); - - if (m_verbosity_level > 0) { - SW_LOG("Finished registration based initialization"); - } -} +void Optimize::InitializeFromRegistration() { RegistrationInitializer(*this).Run(); } //--------------------------------------------------------------------------- void Optimize::RunOptimize() { @@ -2356,7 +1955,8 @@ void Optimize::SetIterationCallback() { this->m_iteration_count++; for (int d = 0; d < m_domains_per_shape; d++) { - current_particle_iterations_ += m_sampler->GetParticleSystem()->GetNumberOfParticles(m_progress_domain_offset + d); + current_particle_iterations_ += + m_sampler->GetParticleSystem()->GetNumberOfParticles(m_progress_domain_offset + d); } if (this->GetShowVisualizer()) { @@ -2614,7 +2214,7 @@ void Optimize::ComputeTotalIterations() { if (m_initialization_mode == InitializationMode::Registration) { const int transfers = std::max(0, GetNumberOfSubjects() - 1) * static_cast(m_domains_per_shape); if (transfers > 0) { - const double share = kTransferProgressShare / (1.0 - kTransferProgressShare); + const double share = TRANSFER_PROGRESS_SHARE / (1.0 - TRANSFER_PROGRESS_SHARE); const auto budget = static_cast(total_particle_iterations_ * share); m_transfer_iteration_weight = static_cast(budget / transfers); total_particle_iterations_ += m_transfer_iteration_weight * transfers; diff --git a/Libs/Optimize/Optimize.h b/Libs/Optimize/Optimize.h index 298aff86efe..b64cf536027 100644 --- a/Libs/Optimize/Optimize.h +++ b/Libs/Optimize/Optimize.h @@ -28,6 +28,8 @@ namespace shapeworks { +class RegistrationInitializer; + class Project; class ParticleGoodBadAssessment; @@ -394,40 +396,17 @@ class Optimize { void RunOptimize(); //! Establish the initial correspondence by spreading particles over a single reference shape and - //! then registering that shape to each of the others + //! then registering that shape to each of the others. RegistrationInitializer does the work; it + //! reaches into this class because the particles, the domains and the settings it needs all live + //! here, and moving them would be a larger change than the one it is worth. void InitializeFromRegistration(); - - //! Pick the shape used as the registration template, honoring an explicitly requested one - int ResolveRegistrationReference(); - - //! Split and optimize particles on the reference shape alone, until it holds the requested counts - void SpreadParticlesOnReference(int reference_shape); - - //! Register the reference to every other shape and carry its particles across - void TransferParticlesFromReference(int reference_shape); + friend class RegistrationInitializer; //! Give the host a chance to refresh during the transfer phase, which runs no optimizer iterations //! and so would otherwise not drive any of the normal per-iteration UI updates. The base //! implementation does nothing; Studio overrides it to redraw the particles and status. virtual void RefreshDuringTransfer() {} - //! Log how well a shape's transferred particles landed on its surface - void ReportTransferQuality(int domain, const std::vector& reference_points, - const std::vector& transferred); - - //! Path a registration's transform is cached at, or "" when caching is off. The name covers - //! everything the transform depends on, so a stale one is never mistaken for a hit. - std::string GetRegistrationCachePath(int reference_domain, int domain) const; - - //! Build the image used to register the given domain. Mesh domains are rasterized to a distance - //! transform; image domains are read back from their groomed path. - Image GetRegistrationImage(int domain); - - //! Return the band to retain around the surface, defaulting to a few voxels of the given spacing - double GetRegistrationBand(double spacing) const; - - //! Return the surface of the given domain, reading it back from disk for image domains - Mesh GetDomainSurface(int domain); //! Return the number of shapes (subjects), as opposed to the total number of domains int GetNumberOfSubjects() const; diff --git a/Libs/Optimize/RegistrationInitializer.cpp b/Libs/Optimize/RegistrationInitializer.cpp new file mode 100644 index 00000000000..a6940b1e534 --- /dev/null +++ b/Libs/Optimize/RegistrationInitializer.cpp @@ -0,0 +1,596 @@ +#include "RegistrationInitializer.h" + +#include +#include +#include + +#include + +#include +#include +#include + +#include "Libs/Optimize/Domain/MeshDomain.h" +#include "Optimize.h" +#include "ShapeworksUtils.h" + +namespace shapeworks { + +//! Width of the retained band around the surface, as a fraction of the shape's largest dimension, +//! when it is not set explicitly. The band is what the similarity metric can actually see, so it +//! must be a physical width independent of the rasterization grid: tying it to the grid would +//! silently narrow it as the grid is raised and starve the metric of context on dissimilar shapes. +//! ~5% of the shape size was the best compromise across similar and dissimilar pairs in testing -- +//! wide enough to bridge shape differences, not so wide that it loses surface detail. +constexpr double REGISTRATION_BAND_FRACTION = 0.05; + +//! Width of the image-domain band, in voxels of the (fixed) distance transform. Image inputs already +//! carry their own resolution, so there is no grid to decouple from. +constexpr double REGISTRATION_BAND_VOXELS = 4.0; + +//! A transferred particle further than this fraction of the shape's size from the surface is +//! considered mislanded, whatever units the data is in +constexpr double TRANSFER_FAR_FRACTION = 0.05; + +//! Warn that a registration may have failed once this fraction of a shape's particles are mislanded. +//! Keyed off a fraction rather than the single worst particle so a stray outlier is not a false alarm. +constexpr double TRANSFER_FAR_FRACTION_THRESHOLD = 0.25; + +//--------------------------------------------------------------------------- +Mesh RegistrationInitializer::GetDomainSurface(int domain) { + auto* particle_domain = optimize_.m_sampler->GetParticleSystem()->GetDomain(domain); + + if (auto* mesh_domain = dynamic_cast(particle_domain)) { + return Mesh(mesh_domain->get_surface()->get_polydata()); + } + + // an image domain keeps only a narrow band of its distance transform, so read the groomed input + // back from disk to recover the full surface + if (domain >= static_cast(optimize_.m_domain_paths.size()) || optimize_.m_domain_paths[domain].empty()) { + throw std::runtime_error( + "Registration based initialization needs the path of each groomed input, but none was given for domain " + + std::to_string(domain)); + } + return Image(optimize_.m_domain_paths[domain]).toMesh(0.0); +} + +//--------------------------------------------------------------------------- +Image RegistrationInitializer::GetRegistrationImage(int domain) { + auto* particle_domain = optimize_.m_sampler->GetParticleSystem()->GetDomain(domain); + + if (particle_domain->GetDomainType() == DomainType::Image) { + if (domain >= static_cast(optimize_.m_domain_paths.size()) || optimize_.m_domain_paths[domain].empty()) { + throw std::runtime_error( + "Registration based initialization needs the path of each groomed input, but none was given for domain " + + std::to_string(domain)); + } + // the groomed input is already a distance transform, so take its resolution as given + Image dt(optimize_.m_domain_paths[domain]); + const auto spacing = dt.spacing(); + const double largest_spacing = std::max({spacing[0], spacing[1], spacing[2]}); + return ImageRegistration::make_registration_image(dt, GetRegistrationBand(largest_spacing)); + } + + // A mesh carries no resolution of its own, so pick one from how large it actually is. Assuming a + // spacing (and a band) in millimeters would break on data stored in other units, where a fixed + // spacing yields either an unusable grid or one too coarse to register. + Mesh mesh = GetDomainSurface(domain); + auto region = mesh.boundingBox(); + const auto extent = region.size(); + const double largest_extent = std::max({extent[0], extent[1], extent[2]}); + const double spacing = largest_extent / optimize_.m_registration_grid_size; + + // Band is a physical width derived from the shape's size, independent of the chosen grid, so the + // grid controls only resolution and the two do not confound each other. + const double band = + optimize_.m_registration_band > 0.0 ? optimize_.m_registration_band : REGISTRATION_BAND_FRACTION * largest_extent; + + // pad far enough that the whole band around the surface stays inside the grid + region.pad(band * 2.0); + return ImageRegistration::make_registration_image( + mesh.toDistanceTransform(region, Point3({spacing, spacing, spacing})), band); +} + +//--------------------------------------------------------------------------- +double RegistrationInitializer::GetRegistrationBand(double spacing) const { + if (optimize_.m_registration_band > 0.0) { + return optimize_.m_registration_band; + } + // default to a band a few voxels wide, which keeps the metric focused near the surface without + // making it so thin that it spans less than a voxel + return REGISTRATION_BAND_VOXELS * spacing; +} + +//--------------------------------------------------------------------------- +int RegistrationInitializer::ResolveRegistrationReference() { + const int num_subjects = optimize_.GetNumberOfSubjects(); + if (num_subjects < 1) { + throw std::runtime_error("No shapes to initialize"); + } + + if (optimize_.m_registration_reference >= 0) { + if (optimize_.m_registration_reference >= num_subjects) { + throw std::runtime_error("Requested registration reference " + + std::to_string(optimize_.m_registration_reference) + + " is out of range, there are only " + std::to_string(num_subjects) + " shapes"); + } + return optimize_.m_registration_reference; + } + + if (num_subjects == 1) { + return 0; + } + + // combine each subject's domains into one mesh so that the template is representative of the whole + // anatomy rather than of a single domain + optimize_.PrintStartMessage("Choosing a registration reference..."); + std::vector meshes; + for (int s = 0; s < num_subjects; s++) { + Mesh mesh = GetDomainSurface(s * optimize_.m_domains_per_shape); + for (int d = 1; d < static_cast(optimize_.m_domains_per_shape); d++) { + mesh += GetDomainSurface(s * optimize_.m_domains_per_shape + d); + } + meshes.push_back(mesh); + } + + const int reference = MeshUtils::findReferenceMesh(meshes); + if (reference < 0 || reference >= num_subjects) { + throw std::runtime_error("Could not choose a registration reference"); + } + optimize_.PrintDoneMessage(); + return reference; +} + +//--------------------------------------------------------------------------- +void RegistrationInitializer::SpreadParticlesOnReference(int reference_shape) { + auto* system = optimize_.m_sampler->GetParticleSystem(); + const int first_domain = reference_shape * optimize_.m_domains_per_shape; + + // only the reference carries particles for the moment, so the correspondence matrices cannot keep + // a consistent layout; they are brought back once every shape has been populated + optimize_.m_sampler->SetCorrespondenceMatricesSuspended(true); + + // seed a single particle on each of the reference's domains + for (int d = 0; d < static_cast(optimize_.m_domains_per_shape); d++) { + const int domain = first_domain + d; + if (system->GetNumberOfParticles(domain) == 0) { + auto* particle_domain = system->GetDomain(domain); + system->AddPosition(particle_domain->GetValidLocationNear(particle_domain->GetZeroCrossingPoint()), domain); + } + } + system->SynchronizePositions(); + + // allocate everything without optimizing, so that turning correspondence off below sticks; Execute + // forces it back on the first time it runs + optimize_.m_sampler->Initialize(); + + // only one shape carries particles at this point, so there is no correspondence to establish yet + optimize_.m_sampler->SetCorrespondenceOff(); + + // report progress against the reference rather than the first shape, which is still empty + optimize_.m_progress_domain_offset = first_domain; + + const double epsilon = optimize_.m_spacing; + + auto needs_split = [&]() { + for (int d = 0; d < static_cast(optimize_.m_domains_per_shape); d++) { + if (system->GetNumberOfParticles(first_domain + d) < optimize_.m_number_of_particles[d]) { + return true; + } + } + return false; + }; + + while (needs_split() && !optimize_.m_aborted) { + optimize_.OptimizerStop(); + + for (int d = 0; d < static_cast(optimize_.m_domains_per_shape); d++) { + const int domain = first_domain + d; + if (system->GetNumberOfParticles(domain) < optimize_.m_number_of_particles[d]) { + system->SplitAllParticlesInDomain(epsilon, domain); + } + } + system->SynchronizePositions(); + + optimize_.m_split_number++; + if (optimize_.m_verbosity_level > 0) { + std::string counts; + for (int d = 0; d < static_cast(optimize_.m_domains_per_shape); d++) { + counts += " " + std::to_string(system->GetNumberOfParticles(first_domain + d)); + } + SW_LOG("Reference split {}, particle count:{}", optimize_.m_split_number, counts); + } + + optimize_.m_energy_a.clear(); + optimize_.m_energy_b.clear(); + optimize_.m_total_energy.clear(); + optimize_.m_str_energy = "split" + std::to_string(optimize_.m_split_number) + "pts_init"; + + optimize_.m_sampler->GetOptimizer()->set_maximum_number_of_iterations(optimize_.m_iterations_per_split); + optimize_.m_sampler->GetOptimizer()->set_number_of_iterations(0); + optimize_.m_sampler->Execute(); + } + + optimize_.m_progress_domain_offset = 0; +} + +//--------------------------------------------------------------------------- +std::string RegistrationInitializer::GetRegistrationCachePath(const ImageRegistration& registration, + int reference_domain, int domain) const { + if (optimize_.m_registration_cache_dir.empty()) { + return ""; + } + + auto describe = [&](int d) { + std::string description = + d < static_cast(optimize_.m_domain_paths.size()) ? optimize_.m_domain_paths[d] : std::to_string(d); + // a groomed input that has been rewritten must not read back as a hit + if (ShapeWorksUtils::file_exists(description)) { + try { + description += ":" + std::to_string(boost::filesystem::file_size(description)); + description += ":" + std::to_string(boost::filesystem::last_write_time(description)); + } catch (const std::exception&) { + } + } + return description; + }; + + // Everything the transform depends on, and nothing that it does not. How the registration itself + // was performed is the registration's own account of it, so that changing the algorithm cannot + // leave a description here saying otherwise; what is left is how these two images were prepared + // for it, which is this class's doing. + std::string key = describe(reference_domain) + "->" + describe(domain); + key += "|" + registration.settings_description(); + key += "|band=" + std::to_string(optimize_.m_registration_band); + key += "|grid=" + std::to_string(optimize_.m_registration_grid_size); + + const auto hash = std::hash{}(key); + std::stringstream name; + name << optimize_.m_registration_cache_dir << "/registration_" << std::hex << hash << ".tfm"; + return name.str(); +} + +//--------------------------------------------------------------------------- +void RegistrationInitializer::TransferParticlesFromReference(int reference_shape) { + auto* system = optimize_.m_sampler->GetParticleSystem(); + const int num_subjects = optimize_.GetNumberOfSubjects(); + + // how each shape ended up landing, so that the run can close by saying what is normal for this + // cohort and which shape sits furthest from it + std::vector> landings; + + for (int d = 0; d < static_cast(optimize_.m_domains_per_shape) && !optimize_.m_aborted; d++) { + const int reference_domain = reference_shape * optimize_.m_domains_per_shape + d; + + std::vector reference_points; + for (auto k = 0; k < system->GetNumberOfParticles(reference_domain); k++) { + reference_points.push_back(system->GetPosition(k, reference_domain)); + } + + // The reference image is only needed to run a registration, so it is built lazily on the first + // cache miss. Rasterizing a large mesh takes several seconds, and when every transfer is a cache + // hit (e.g. re-running with different optimization parameters) it is never needed at all. + std::optional reference_image; + + for (int s = 0; s < num_subjects && !optimize_.m_aborted; s++) { + if (s == reference_shape) { + continue; + } + const int domain = s * optimize_.m_domains_per_shape + d; + + optimize_.UpdateProgress(fmt::format("Registering shape {} of {}", s + 1, num_subjects)); + optimize_.RefreshDuringTransfer(); + + ImageRegistration registration; + registration.set_transform_type(optimize_.m_registration_transform_type); + registration.set_gradient_step(optimize_.m_registration_gradient_step); + registration.set_update_field_variance(optimize_.m_registration_flow_sigma); + + const auto cache_path = GetRegistrationCachePath(registration, reference_domain, domain); + // check for the file first: asking the reader for one that is not there is a normal cache + // miss, but it makes the HDF5 layer print an alarming block of text + const bool cached = + !cache_path.empty() && ShapeWorksUtils::file_exists(cache_path) && registration.load_transform(cache_path); + if (cached) { + SW_DEBUG("Reusing cached registration: {}", cache_path); + } else { + if (!reference_image) { + // only now, on a genuine miss, is the reference image worth building + optimize_.UpdateProgress(optimize_.m_domains_per_shape > 1 + ? fmt::format("Preparing registration (domain {})", d + 1) + : std::string("Preparing registration")); + optimize_.RefreshDuringTransfer(); + reference_image = GetRegistrationImage(reference_domain); + } + registration.run(*reference_image, GetRegistrationImage(domain)); + if (!cache_path.empty()) { + try { + registration.save_transform(cache_path); + } catch (const std::exception& e) { + SW_WARN("Unable to cache registration: {}", e.what()); + } + } + } + + auto transferred = registration.transform_points(reference_points); + + // A registration can leave a particle somewhere this shape cannot be sampled: outside its + // image, or too far from its surface for the narrow band to reach. Settle that here, where + // the shape can be named, rather than leaving it to surface as an unplaceable particle from + // deep inside the sampler. + RescueTransferredParticles(domain, reference_points, transferred); + + // AddPositionList applies the domain constraints, which pulls each point onto the surface + system->AddPositionList(transferred, domain); + + landings.emplace_back(GetDomainName(domain), ReportTransferQuality(domain, reference_points, transferred)); + + // registrations do not run particle iterations, so charge each one its share of the budget + // reserved for them; without this the bar would sit still through the longest phase + optimize_.current_particle_iterations_ += optimize_.m_transfer_iteration_weight; + optimize_.UpdateProgress(fmt::format("Registered shape {} of {}", s + 1, num_subjects)); + // the shape now has its particles; let the GUI redraw them (and the status) between + // registrations rather than only when the whole phase ends + optimize_.RefreshDuringTransfer(); + } + } + + ReportTransferSummary(landings); +} + +//--------------------------------------------------------------------------- +std::string RegistrationInitializer::GetDomainName(int domain) const { + return domain < static_cast(optimize_.m_filenames.size()) ? optimize_.m_filenames[domain] + : std::to_string(domain); +} + +//--------------------------------------------------------------------------- +void RegistrationInitializer::ReportTransferSummary(const std::vector>& landings) const { + if (landings.empty()) { + return; + } + + std::vector distances; + distances.reserve(landings.size()); + for (const auto& landing : landings) { + distances.push_back(landing.second); + } + const auto middle = distances.begin() + distances.size() / 2; + std::nth_element(distances.begin(), middle, distances.end()); + + const auto furthest = std::max_element(landings.begin(), landings.end(), + [](const auto& a, const auto& b) { return a.second < b.second; }); + + // One shape landing further out than the rest is what a failed registration looks like from here, + // so give the typical figure and the outlier together: either number alone says nothing. + SW_LOG( + "Registration based initialization: across {} registered shapes, particles landed a median of {:.3g} from the " + "surface, and furthest on {} at {:.3g}", + landings.size(), *middle, furthest->first, furthest->second); +} + +//--------------------------------------------------------------------------- +void RegistrationInitializer::RescueTransferredParticles(int domain, const std::vector& reference_points, + std::vector& transferred) const { + auto* particle_domain = optimize_.m_sampler->GetParticleSystem()->GetDomain(domain); + + const auto& lower = particle_domain->GetLowerBound(); + const auto& upper = particle_domain->GetUpperBound(); + + std::vector placed(transferred.size(), true); + std::vector misplaced; + int off_the_image = 0; + for (size_t i = 0; i < transferred.size(); i++) { + placed[i] = particle_domain->IsValidLocation(transferred[i]); + if (placed[i]) { + continue; + } + misplaced.push_back(i); + + // Whether a particle can be sampled depends on how wide a narrow band was asked for, which is a + // setting rather than a verdict on the registration; whether it is inside the image at all does + // not. Count the two apart so that a tighter band cannot masquerade as a failed registration. + for (unsigned int c = 0; c < 3; c++) { + if (transferred[i][c] < lower[c] || transferred[i][c] > upper[c]) { + off_the_image++; + break; + } + } + } + + if (misplaced.empty()) { + return; + } + + const std::string name = GetDomainName(domain); + + // A registration that worked puts every particle somewhere on the shape, so more than a stray + // handful landing outside the image altogether means it did not work at all. Nothing downstream + // can make a shape model out of that, and stopping here is what lets the failure name the shape it + // happened to. + if (off_the_image > TRANSFER_FAR_FRACTION_THRESHOLD * transferred.size()) { + throw std::runtime_error(fmt::format( + "Registration based initialization failed for {}: {} of its {} particles were mapped outside its image " + "altogether. The registration from the reference shape did not converge, so there is nothing here to start " + "the optimization from. Check that this shape's groomed input covers the same anatomy as the rest, or " + "choose a different reference shape with the initialization reference setting.", + name, off_the_image, transferred.size())); + } + + // Aim each stray at where its neighbours ended up rather than at some arbitrary point on the + // shape. Particles are transferred in correspondence, so the one beside it on the reference is + // the best evidence there is of where this one belonged. + constexpr int RESCUE_STEPS = 40; + + double furthest = 0.0; + + for (const auto index : misplaced) { + // the nearest particle that was placed, measured on the reference, where the correspondence + // between the two shapes is defined + size_t nearest = 0; + double nearest_distance = std::numeric_limits::max(); + for (size_t i = 0; i < transferred.size(); i++) { + const double distance = reference_points[index].SquaredEuclideanDistanceTo(reference_points[i]); + if (placed[i] && distance < nearest_distance) { + nearest_distance = distance; + nearest = i; + } + } + + const Point3 stray = transferred[index]; + const Point3 target = transferred[nearest]; + auto along = [&](double fraction) { + Point3 point; + for (unsigned int c = 0; c < 3; c++) { + point[c] = stray[c] + fraction * (target[c] - stray[c]); + } + return point; + }; + + // Walk in from the stray until the shape can be sampled again, then carry on halfway to the + // target and settle there. The first location that can be sampled is the far edge of the narrow + // band, where there is too little distance transform left around the particle to project it onto + // the surface; going halfway on towards a particle that did land keeps the direction the + // registration chose while leaving this one somewhere the domain can work with. The last step + // of the walk is the target itself, so somewhere is always found. + for (int step = 1; step <= RESCUE_STEPS; step++) { + const double fraction = static_cast(step) / RESCUE_STEPS; + if (!particle_domain->IsValidLocation(along(fraction))) { + continue; + } + const Point3 deeper = along((fraction + 1.0) / 2.0); + transferred[index] = particle_domain->IsValidLocation(deeper) ? deeper : along(fraction); + break; + } + + furthest = std::max(furthest, stray.EuclideanDistanceTo(transferred[index])); + } + + // Telling someone that something is worth checking is no use without telling them how. Which + // remedy to point at depends on where the particles went: outside the image is the registration's + // doing, whereas inside it but beyond the narrow band only means the band is too tight to hold + // them, which says nothing about the registration and has an entirely different answer. + if (off_the_image > 0) { + SW_WARN( + "{}: {} of {} transferred particles landed off the shape, by as much as {:.3g}, and were moved back onto it. " + "The registration from the reference shape did not fit this one closely (the distance a well registered shape " + "lands within is logged at the end of initialization). Look at this shape in the viewer, or in the " + "correspondence quality panel, once the run finishes. If its particles do not correspond to the others, " + "either register from a different reference shape with the initialization reference setting, or check that " + "grooming left this shape aligned with the rest of the cohort.", + name, misplaced.size(), transferred.size(), furthest); + } else { + SW_WARN( + "{}: {} of {} transferred particles registered within the image but outside the narrow band, by as much as " + "{:.3g}, and were moved back onto the shape. Increasing the narrow band optimization parameter may help.", + name, misplaced.size(), transferred.size(), furthest); + } +} + +//--------------------------------------------------------------------------- +double RegistrationInitializer::ReportTransferQuality(int domain, const std::vector& reference_points, + const std::vector& transferred) { + if (transferred.empty()) { + return 0.0; + } + + auto* system = optimize_.m_sampler->GetParticleSystem(); + + // Judge the transfer against the size of the shape rather than an absolute distance, so the same + // thresholds work whatever units and resolution the data is in. Take that size from the domain + // itself, not from the particles: particles that have collapsed together would otherwise shrink + // the scale in step with the distances being judged, and always look acceptable. + const auto* particle_domain = system->GetDomain(domain); + const double shape_scale = particle_domain->GetLowerBound().EuclideanDistanceTo(particle_domain->GetUpperBound()); + + // A particle beyond this fraction of the shape's size is clearly in the wrong place + const double far_distance = TRANSFER_FAR_FRACTION * shape_scale; + + double total_snap = 0.0; + double worst_snap = 0.0; + int far_count = 0; + int unmoved = 0; + + for (size_t i = 0; i < transferred.size(); i++) { + // how far the point had to travel to reach the surface is how far off the surface it landed + const double snap = transferred[i].EuclideanDistanceTo(system->GetPosition(i, domain)); + total_snap += snap; + worst_snap = std::max(worst_snap, snap); + if (shape_scale > 0.0 && snap > far_distance) { + far_count++; + } + + // a point outside the displacement field is returned unchanged rather than reported as an error, + // so an unmoved point means the field did not cover it + if (transferred[i] == reference_points[i]) { + unmoved++; + } + } + + const double mean_snap = total_snap / transferred.size(); + const double far_fraction = static_cast(far_count) / transferred.size(); + // a registration that failed leaves many particles far from the surface, not just one outlier + const bool suspect = far_fraction > TRANSFER_FAR_FRACTION_THRESHOLD || unmoved > 0; + + const std::string name = GetDomainName(domain); + + if (optimize_.m_verbosity_level > 0 || suspect) { + SW_LOG("{}: transferred particles landed {:.3g} from the surface on average (worst {:.3g})", name, mean_snap, + worst_snap); + } + + if (unmoved > 0) { + SW_WARN("{}: {} of {} transferred particles fell outside the registration field and did not move", name, unmoved, + transferred.size()); + } + + if (far_fraction > TRANSFER_FAR_FRACTION_THRESHOLD) { + SW_WARN("{}: {:.0f}% of transferred particles landed far from the surface, the registration may have failed", name, + 100.0 * far_fraction); + } + + return mean_snap; +} + +//--------------------------------------------------------------------------- +void RegistrationInitializer::Run() { + const int reference = ResolveRegistrationReference(); + optimize_.m_registration_reference_chosen = reference; + + const int first_domain = reference * optimize_.m_domains_per_shape; + const std::string name = first_domain < static_cast(optimize_.m_filenames.size()) + ? optimize_.m_filenames[first_domain] + : std::to_string(reference); + SW_LOG("Spreading particles on reference shape {} ({})", reference, name); + + SpreadParticlesOnReference(reference); + + // each shape is populated with a full set below, so the matrices can track them again + optimize_.m_sampler->SetCorrespondenceMatricesSuspended(false); + + if (!optimize_.m_aborted) { + TransferParticlesFromReference(reference); + } + + // Only now does every shape hold its particles. The matrices size themselves from the first + // shape's domains, so they cannot be brought up to date any earlier: while particles were being + // spread, only the reference (which is usually not the first shape) had any. + auto* system = optimize_.m_sampler->GetParticleSystem(); + system->ResyncObservers(); + system->SynchronizePositions(); + + // every shape now carries a full set of corresponding particles + optimize_.m_sampler->SetCorrespondenceOn(); + + optimize_.WritePointFiles(); + optimize_.WritePointFilesWithFeatures(); + optimize_.WriteTransformFile(); + optimize_.WriteTransformFiles(); + optimize_.WriteCuttingPlanePoints(); + + if (optimize_.m_verbosity_level > 0) { + SW_LOG("Finished registration based initialization"); + } +} + +} // namespace shapeworks diff --git a/Libs/Optimize/RegistrationInitializer.h b/Libs/Optimize/RegistrationInitializer.h new file mode 100644 index 00000000000..1b5ae5c1500 --- /dev/null +++ b/Libs/Optimize/RegistrationInitializer.h @@ -0,0 +1,81 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +#include "Libs/Common/Shapeworks.h" + +namespace shapeworks { + +class Optimize; + +/** + * \class RegistrationInitializer + * \ingroup Group-Optimize + * + * Establishes the initial correspondence by a different route than splitting. Particles are spread + * over one reference shape alone, that shape is registered to each of the others, and its particles + * are carried across, so that every shape starts from a set that already corresponds rather than + * from one that has to be brought into correspondence by optimization. + * + * Holds no state of its own beyond the run in progress: the particles, the domains and the settings + * all belong to the Optimize it is initializing, which it works on directly. + */ +class RegistrationInitializer { + public: + explicit RegistrationInitializer(Optimize& optimize) : optimize_(optimize) {} + + //! Spread particles over the reference shape and carry them to every other shape + void Run(); + + private: + //! Pick the shape used as the registration template, honoring an explicitly requested one + int ResolveRegistrationReference(); + + //! Split and optimize particles on the reference shape alone, until it holds the requested counts + void SpreadParticlesOnReference(int reference_shape); + + //! Register the reference to every other shape and carry its particles across + void TransferParticlesFromReference(int reference_shape); + + //! Move transferred particles that the domain cannot sample back to somewhere it can, and refuse + //! a registration that has misplaced more than a stray handful of them + void RescueTransferredParticles(int domain, const std::vector& reference_points, + std::vector& transferred) const; + + //! Log how well a shape's transferred particles landed on its surface, and return how far from it + //! they landed on average, so that the shapes can be compared with one another + double ReportTransferQuality(int domain, const std::vector& reference_points, + const std::vector& transferred); + + //! Log how the whole cohort landed, which is what gives any one shape's figure a scale to be read + //! against + void ReportTransferSummary(const std::vector>& landings) const; + + //! The groomed file a domain came from, for messages that need to name it + std::string GetDomainName(int domain) const; + + //! Path the given registration's transform is cached at, or "" when caching is off. The name + //! covers everything the transform depends on, so a stale one is never mistaken for a hit. + std::string GetRegistrationCachePath(const ImageRegistration& registration, int reference_domain, + int domain) const; + + //! Build the image used to register the given domain. Mesh domains are rasterized to a distance + //! transform; image domains are read back from their groomed path. + Image GetRegistrationImage(int domain); + + //! Return the band to retain around the surface, defaulting to a few voxels of the given spacing + double GetRegistrationBand(double spacing) const; + + //! Return the surface of the given domain, reading it back from disk for image domains + Mesh GetDomainSurface(int domain); + + Optimize& optimize_; +}; + +} // namespace shapeworks From 3c6868b98e29af0ac66f7d9eaaf722e7cc3921fe Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Wed, 26 Aug 2026 12:34:14 -0600 Subject: [PATCH 2/5] Allow copying messages out of the Studio message history The list was built with selection and keyboard focus disabled, so the only way to get the text of an error out of Studio was to retype it from the screen. Messages can now be selected and copied with the usual shortcut or from a context menu, and an arriving message scrolls into view without clearing a selection that is part way through being made. --- Studio/Interface/LogWindow.cpp | 50 +++++++++++++++++++++++++++++++++- Studio/Interface/LogWindow.h | 12 ++++++++ Studio/Interface/LogWindow.ui | 9 ++++-- 3 files changed, 67 insertions(+), 4 deletions(-) diff --git a/Studio/Interface/LogWindow.cpp b/Studio/Interface/LogWindow.cpp index f4f498fcec0..547120a2bb0 100644 --- a/Studio/Interface/LogWindow.cpp +++ b/Studio/Interface/LogWindow.cpp @@ -1,5 +1,10 @@ #include "LogWindow.h" +#include +#include +#include +#include + // Automatically generated UI file #include "ui_LogWindow.h" @@ -11,6 +16,19 @@ LogWindow::LogWindow(QWidget* parent) : QDialog(parent) { ui_->setupUi(this); ui_->history_list_->setWordWrap(true); + // The message history is where someone goes to get the text of an error back out of Studio, so it + // has to be possible to take it: select the lines that matter, or copy the lot. + auto* copy = new QAction(tr("&Copy"), ui_->history_list_); + copy->setShortcut(QKeySequence::Copy); + // window rather than widget scope, so the shortcut works without first clicking into the list + copy->setShortcutContext(Qt::WindowShortcut); + connect(copy, &QAction::triggered, this, &LogWindow::copy_selected_messages); + ui_->history_list_->addAction(copy); + + auto* copy_all = new QAction(tr("Copy &All"), ui_->history_list_); + connect(copy_all, &QAction::triggered, this, &LogWindow::copy_all_messages); + ui_->history_list_->addAction(copy_all); + QIcon icon = windowIcon(); Qt::WindowFlags flags = windowFlags(); Qt::WindowFlags helpFlag = Qt::WindowContextHelpButtonHint; @@ -43,7 +61,37 @@ void LogWindow::add_message(MessageType message_type, QString message) { QListWidgetItem* item = new QListWidgetItem(message, this->ui_->history_list_); item->setForeground(color); this->ui_->history_list_->addItem(item); - this->ui_->history_list_->setCurrentItem(item); + // follow the newest message without touching the selection, which someone may be part way through + // making when the next message arrives + this->ui_->history_list_->scrollToItem(item); +} + +//--------------------------------------------------------------------------- +void LogWindow::copy_selected_messages() { + auto items = ui_->history_list_->selectedItems(); + if (items.isEmpty()) { + copy_all_messages(); + return; + } + copy_messages(items); +} + +//--------------------------------------------------------------------------- +void LogWindow::copy_all_messages() { + QList items; + for (int i = 0; i < ui_->history_list_->count(); i++) { + items.append(ui_->history_list_->item(i)); + } + copy_messages(items); +} + +//--------------------------------------------------------------------------- +void LogWindow::copy_messages(const QList& items) { + QStringList messages; + for (const auto* item : items) { + messages << item->text(); + } + QApplication::clipboard()->setText(messages.join("\n")); } } // namespace shapeworks diff --git a/Studio/Interface/LogWindow.h b/Studio/Interface/LogWindow.h index 0fc008051a3..f452e20c9a7 100644 --- a/Studio/Interface/LogWindow.h +++ b/Studio/Interface/LogWindow.h @@ -1,6 +1,7 @@ #pragma once #include +#include #include @@ -24,8 +25,19 @@ Q_OBJECT //! add a message to the history void add_message(MessageType message_type, QString message); +private Q_SLOTS: + + //! put the selected messages on the clipboard, or every message when none are selected + void copy_selected_messages(); + + //! put every message on the clipboard + void copy_all_messages(); + private: + //! put the given messages on the clipboard, one per line + void copy_messages(const QList& items); + Ui_LogWindow* ui_; }; diff --git a/Studio/Interface/LogWindow.ui b/Studio/Interface/LogWindow.ui index 6f548f8011e..4015ffbde27 100644 --- a/Studio/Interface/LogWindow.ui +++ b/Studio/Interface/LogWindow.ui @@ -83,10 +83,10 @@ QLabel#warning_label_{ - Qt::NoFocus + Qt::StrongFocus - Qt::DefaultContextMenu + Qt::ActionsContextMenu 1 @@ -104,7 +104,10 @@ QLabel#warning_label_{ true - QAbstractItemView::NoSelection + QAbstractItemView::ExtendedSelection + + + Qt::ElideNone QAbstractItemView::ScrollPerItem From b3b86498ce96cc1461ec2b825f9a2a603e4e675e Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Wed, 26 Aug 2026 13:00:26 -0600 Subject: [PATCH 3/5] Address review of the registration initialization changes The rescue walks a stray particle towards one that landed, so a shape with no placed particle at all left every stray where it was and then threw from inside the sampler, which is the failure the rescue exists to prevent; that state now reports itself as a band too narrow to hold the registration. ITK throws rather than returning a zero mass for an empty image, so the scale seed's guard against one was unreachable and cost the pair its affine stage. A transform file is written after its displacement field rather than before, so that a run interrupted between the two leaves no entry instead of one that reads back as a deformable registration quietly reduced to a linear one. The narrow band is quoted in the voxels the parameter is set in rather than the world units the domain keeps it in, the affine seed no longer sets a centre that ImageRegistrationMethodv4 overrides, and the transfer summary counts registrations rather than shapes now that it can hold more than one per shape. --- Libs/Image/ImageRegistration.cpp | 44 ++++++++++++++--------- Libs/Optimize/Domain/ImageDomain.h | 12 ++++--- Libs/Optimize/Optimize.cpp | 3 -- Libs/Optimize/RegistrationInitializer.cpp | 33 +++++++++++------ 4 files changed, 58 insertions(+), 34 deletions(-) diff --git a/Libs/Image/ImageRegistration.cpp b/Libs/Image/ImageRegistration.cpp index ef5716edf91..5a058110604 100644 --- a/Libs/Image/ImageRegistration.cpp +++ b/Libs/Image/ImageRegistration.cpp @@ -105,12 +105,16 @@ std::vector clamp_shrink_factors(const std::vector& struct ScaleSeed { bool usable{false}; double scale{1.0}; - ImageType::PointType center; }; ScaleSeed estimate_scale_seed(const ImageType* fixed, const ImageType* moving) { using MomentsType = itk::ImageMomentsCalculator; + ScaleSeed seed; + + // An image with nothing in it has no moments to take, and ITK says so by throwing rather than by + // returning a mass of zero. There is no size to compare in that case, which is a reason to leave + // the stage as it was and not a reason to fail the registration. auto measure = [](const ImageType* image) { auto moments = MomentsType::New(); moments->SetImage(const_cast(image)); @@ -118,9 +122,14 @@ ScaleSeed estimate_scale_seed(const ImageType* fixed, const ImageType* moving) { return moments; }; - ScaleSeed seed; - auto fixed_moments = measure(fixed); - auto moving_moments = measure(moving); + MomentsType::Pointer fixed_moments; + MomentsType::Pointer moving_moments; + try { + fixed_moments = measure(fixed); + moving_moments = measure(moving); + } catch (const itk::ExceptionObject&) { + return seed; + } const double fixed_mass = fixed_moments->GetTotalMass(); const double moving_mass = moving_moments->GetTotalMass(); @@ -151,10 +160,6 @@ ScaleSeed estimate_scale_seed(const ImageType* fixed, const ImageType* moving) { return seed; } - const auto center_of_gravity = fixed_moments->GetCenterOfGravity(); - for (unsigned int i = 0; i < 3; i++) { - seed.center[i] = center_of_gravity[i]; - } seed.scale = scale; seed.usable = true; return seed; @@ -401,13 +406,14 @@ void ImageRegistration::Impl::run_affine() { auto metric = LinearMetricType::New(); using RegistrationType = itk::ImageRegistrationMethodv4; - // Start from the size difference between the two shapes rather than from the identity. The - // composite applies what it holds in reverse, so this affine acts on points still in fixed space, - // and the point to scale about is the fixed image's own centre of gravity. + // Start from the size difference between the two shapes rather than from the identity. Only the + // matrix is set: ImageRegistrationMethodv4 initializes the centre of a linear output transform + // itself, from the last linear transform of the moving initial transform, so the scaling pivots + // about the rigid stage's centre whatever is set here. That is the fixed image's centre of + // gravity, which is the point to scale about anyway. auto affine = AffineTransformType::New(); const auto seed = estimate_scale_seed(fixed, moving); if (seed.usable) { - affine->SetCenter(seed.center); affine->Scale(seed.scale); } @@ -720,11 +726,10 @@ void ImageRegistration::save_transform(const std::string& filename) const { } try { - auto writer = itk::TransformFileWriterTemplate::New(); - writer->SetFileName(filename); - writer->SetInput(linear); - writer->Update(); - + // The field goes first and the transform file last, because load_transform takes the transform + // file as the record that an entry exists and treats a missing field as a registration that had + // none. Written the other way round, a run interrupted between the two leaves an entry that + // reads back as a deformable registration silently downgraded to a linear one. if (stored_field) { auto field_writer = itk::ImageFileWriter::New(); field_writer->SetFileName(displacement_field_path(filename)); @@ -732,6 +737,11 @@ void ImageRegistration::save_transform(const std::string& filename) const { field_writer->UseCompressionOn(); field_writer->Update(); } + + auto writer = itk::TransformFileWriterTemplate::New(); + writer->SetFileName(filename); + writer->SetInput(linear); + writer->Update(); } catch (const itk::ExceptionObject& e) { throw std::runtime_error(std::string("unable to write transform \"") + filename + "\": " + e.what()); } diff --git a/Libs/Optimize/Domain/ImageDomain.h b/Libs/Optimize/Domain/ImageDomain.h index a6431286c42..9e6d03c64b2 100644 --- a/Libs/Optimize/Domain/ImageDomain.h +++ b/Libs/Optimize/Domain/ImageDomain.h @@ -232,11 +232,11 @@ class ImageDomain : public ParticleRegionDomain { << ". "; // a particle is clamped to the bounding box before it is sampled, so one that arrived from - // outside the image is now sitting exactly on one of its faces - const double tolerance = GetSpacing().GetVnlVector().max_value(); + // outside the image is now sitting within a voxel of one of its faces + const double voxel = GetSpacing().GetVnlVector().max_value(); bool on_edge = false; for (unsigned int i = 0; i < DIMENSION; i++) { - on_edge = on_edge || p[i] <= GetLowerBound()[i] + tolerance || p[i] >= GetUpperBound()[i] - tolerance; + on_edge = on_edge || p[i] <= GetLowerBound()[i] + voxel || p[i] >= GetUpperBound()[i] - voxel; } if (on_edge) { @@ -246,7 +246,11 @@ class ImageDomain : public ParticleRegionDomain { "begin with, usually by an initialization that did not converge or by shapes that grooming left " "unaligned."; } else { - message << "It is inside the image but further than the narrow band (" << m_NarrowBand + // quoted in voxels, which is what the narrow band optimization parameter is set in; the domain + // keeps it as a width in world units, and telling someone to raise a number they never typed + // has them lower the one they did + message << "It is inside the image but further than the narrow band (" << m_NarrowBand / voxel + << " voxels, " << m_NarrowBand << ") from the surface, so no distance transform is kept there. Consider increasing the narrow band " "optimization parameter."; } diff --git a/Libs/Optimize/Optimize.cpp b/Libs/Optimize/Optimize.cpp index f28da76e98f..5aa82a1d4ab 100644 --- a/Libs/Optimize/Optimize.cpp +++ b/Libs/Optimize/Optimize.cpp @@ -2,10 +2,7 @@ #include #include #include -#include -#include // #include -#include #include #include #include diff --git a/Libs/Optimize/RegistrationInitializer.cpp b/Libs/Optimize/RegistrationInitializer.cpp index a6940b1e534..9a17bf545a9 100644 --- a/Libs/Optimize/RegistrationInitializer.cpp +++ b/Libs/Optimize/RegistrationInitializer.cpp @@ -2,13 +2,13 @@ #include #include +#include #include #include #include #include -#include #include "Libs/Optimize/Domain/MeshDomain.h" #include "Optimize.h" @@ -363,15 +363,19 @@ void RegistrationInitializer::ReportTransferSummary(const std::vectorfirst, furthest->second); } //--------------------------------------------------------------------------- -void RegistrationInitializer::RescueTransferredParticles(int domain, const std::vector& reference_points, - std::vector& transferred) const { +void RegistrationInitializer::RescueTransferredParticles(int domain, + const std::vector& reference_points, + std::vector& transferred) const { auto* particle_domain = optimize_.m_sampler->GetParticleSystem()->GetDomain(domain); const auto& lower = particle_domain->GetLowerBound(); @@ -417,6 +421,18 @@ void RegistrationInitializer::RescueTransferredParticles(int domain, const std:: name, off_the_image, transferred.size())); } + // The rescue walks each stray towards a particle that did land, so there has to be one. Every + // particle being unplaceable means the band is too narrow to hold this registration at all, which + // no amount of walking will fix, and walking towards a stray would leave them where they were. + if (misplaced.size() == transferred.size()) { + throw std::runtime_error(fmt::format( + "Registration based initialization failed for {}: none of its {} particles landed within the narrow band of " + "the surface, so there is nowhere to move them back to. Registration based initialization starts particles " + "wherever the reference shape maps to, which can be much further from the surface than the split based " + "initialization the default narrow band is sized for; increase the narrow band optimization parameter.", + name, transferred.size())); + } + // Aim each stray at where its neighbours ended up rather than at some arbitrary point on the // shape. Particles are transferred in correspondence, so the one beside it on the reference is // the best evidence there is of where this one belonged. @@ -489,7 +505,7 @@ void RegistrationInitializer::RescueTransferredParticles(int domain, const std:: //--------------------------------------------------------------------------- double RegistrationInitializer::ReportTransferQuality(int domain, const std::vector& reference_points, - const std::vector& transferred) { + const std::vector& transferred) { if (transferred.empty()) { return 0.0; } @@ -557,11 +573,8 @@ void RegistrationInitializer::Run() { const int reference = ResolveRegistrationReference(); optimize_.m_registration_reference_chosen = reference; - const int first_domain = reference * optimize_.m_domains_per_shape; - const std::string name = first_domain < static_cast(optimize_.m_filenames.size()) - ? optimize_.m_filenames[first_domain] - : std::to_string(reference); - SW_LOG("Spreading particles on reference shape {} ({})", reference, name); + SW_LOG("Spreading particles on reference shape {} ({})", reference, + GetDomainName(reference * optimize_.m_domains_per_shape)); SpreadParticlesOnReference(reference); From 39884461b3dec87f3393ad488c6869dd427860a9 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Wed, 26 Aug 2026 13:20:38 -0600 Subject: [PATCH 4/5] Judge a transferred particle by the point the domain will sample A particle is clamped into the domain's bounding box before it is sampled, but it was being judged where the registration left it, so one a fraction of a voxel outside a shape tight against its image counted against the threshold that abandons the run, for a registration that clamping would have made usable. The walk that moves a stray back now asks the same question, since a target that only clamping makes good is one it could otherwise never reach. The check for particles the displacement field did not cover is gone. It tested transferred and reference points for exact equality, which the rigid stage in the composite makes impossible, and the field is defined over the fixed image, which the reference particles are always well inside; so it guarded a state that cannot arise with a test that could not detect it. --- Libs/Optimize/RegistrationInitializer.cpp | 37 +++++++++++------------ Libs/Optimize/RegistrationInitializer.h | 3 +- 2 files changed, 19 insertions(+), 21 deletions(-) diff --git a/Libs/Optimize/RegistrationInitializer.cpp b/Libs/Optimize/RegistrationInitializer.cpp index 9a17bf545a9..4fbb1dd9b78 100644 --- a/Libs/Optimize/RegistrationInitializer.cpp +++ b/Libs/Optimize/RegistrationInitializer.cpp @@ -323,7 +323,7 @@ void RegistrationInitializer::TransferParticlesFromReference(int reference_shape // AddPositionList applies the domain constraints, which pulls each point onto the surface system->AddPositionList(transferred, domain); - landings.emplace_back(GetDomainName(domain), ReportTransferQuality(domain, reference_points, transferred)); + landings.emplace_back(GetDomainName(domain), ReportTransferQuality(domain, transferred)); // registrations do not run particle iterations, so charge each one its share of the budget // reserved for them; without this the bar would sit still through the longest phase @@ -381,11 +381,23 @@ void RegistrationInitializer::RescueTransferredParticles(int domain, const auto& lower = particle_domain->GetLowerBound(); const auto& upper = particle_domain->GetUpperBound(); + // The domain clamps a particle into its bounding box before sampling it, so ask whether the point + // the domain will actually see can be sampled rather than the one the registration produced: a + // particle a fraction of a voxel outside the image is not misplaced at all if clamping puts it + // back within the band. + auto samplable = [&](const Point3& point) { + Point3 clamped = point; + for (unsigned int c = 0; c < 3; c++) { + clamped[c] = std::clamp(clamped[c], lower[c], upper[c]); + } + return particle_domain->IsValidLocation(clamped); + }; + std::vector placed(transferred.size(), true); std::vector misplaced; int off_the_image = 0; for (size_t i = 0; i < transferred.size(); i++) { - placed[i] = particle_domain->IsValidLocation(transferred[i]); + placed[i] = samplable(transferred[i]); if (placed[i]) { continue; } @@ -471,11 +483,11 @@ void RegistrationInitializer::RescueTransferredParticles(int domain, // of the walk is the target itself, so somewhere is always found. for (int step = 1; step <= RESCUE_STEPS; step++) { const double fraction = static_cast(step) / RESCUE_STEPS; - if (!particle_domain->IsValidLocation(along(fraction))) { + if (!samplable(along(fraction))) { continue; } const Point3 deeper = along((fraction + 1.0) / 2.0); - transferred[index] = particle_domain->IsValidLocation(deeper) ? deeper : along(fraction); + transferred[index] = samplable(deeper) ? deeper : along(fraction); break; } @@ -504,8 +516,7 @@ void RegistrationInitializer::RescueTransferredParticles(int domain, } //--------------------------------------------------------------------------- -double RegistrationInitializer::ReportTransferQuality(int domain, const std::vector& reference_points, - const std::vector& transferred) { +double RegistrationInitializer::ReportTransferQuality(int domain, const std::vector& transferred) { if (transferred.empty()) { return 0.0; } @@ -525,7 +536,6 @@ double RegistrationInitializer::ReportTransferQuality(int domain, const std::vec double total_snap = 0.0; double worst_snap = 0.0; int far_count = 0; - int unmoved = 0; for (size_t i = 0; i < transferred.size(); i++) { // how far the point had to travel to reach the surface is how far off the surface it landed @@ -535,18 +545,12 @@ double RegistrationInitializer::ReportTransferQuality(int domain, const std::vec if (shape_scale > 0.0 && snap > far_distance) { far_count++; } - - // a point outside the displacement field is returned unchanged rather than reported as an error, - // so an unmoved point means the field did not cover it - if (transferred[i] == reference_points[i]) { - unmoved++; - } } const double mean_snap = total_snap / transferred.size(); const double far_fraction = static_cast(far_count) / transferred.size(); // a registration that failed leaves many particles far from the surface, not just one outlier - const bool suspect = far_fraction > TRANSFER_FAR_FRACTION_THRESHOLD || unmoved > 0; + const bool suspect = far_fraction > TRANSFER_FAR_FRACTION_THRESHOLD; const std::string name = GetDomainName(domain); @@ -555,11 +559,6 @@ double RegistrationInitializer::ReportTransferQuality(int domain, const std::vec worst_snap); } - if (unmoved > 0) { - SW_WARN("{}: {} of {} transferred particles fell outside the registration field and did not move", name, unmoved, - transferred.size()); - } - if (far_fraction > TRANSFER_FAR_FRACTION_THRESHOLD) { SW_WARN("{}: {:.0f}% of transferred particles landed far from the surface, the registration may have failed", name, 100.0 * far_fraction); diff --git a/Libs/Optimize/RegistrationInitializer.h b/Libs/Optimize/RegistrationInitializer.h index 1b5ae5c1500..494651faf43 100644 --- a/Libs/Optimize/RegistrationInitializer.h +++ b/Libs/Optimize/RegistrationInitializer.h @@ -50,8 +50,7 @@ class RegistrationInitializer { //! Log how well a shape's transferred particles landed on its surface, and return how far from it //! they landed on average, so that the shapes can be compared with one another - double ReportTransferQuality(int domain, const std::vector& reference_points, - const std::vector& transferred); + double ReportTransferQuality(int domain, const std::vector& transferred); //! Log how the whole cohort landed, which is what gives any one shape's figure a scale to be read //! against From 1259ada4735cc472043bc4f0461f67c968d0cb71 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Wed, 26 Aug 2026 14:38:05 -0600 Subject: [PATCH 5/5] Judge the overlap guard on the step that was taken The guard read the number of samples the metric last found usable, but with a line search that count belongs to whichever trial the golden section search evaluated last rather than to the step actually taken, so it could condemn a good step and, worse, pass a bad one. It now maps a few hundred points of its own through the transform and counts how many land on the moving image, which measures the step that was taken and stays comparable across levels of the pyramid rather than growing with them. Two tests cover what the registration gained: that the affine stage recovers a size difference too large for it to find unaided, which fails without the seed that hands it one, and that every setting able to change a transform also changes the description a caller keeps transforms against. --- Libs/Image/ImageRegistration.cpp | 74 ++++++++++++++++++++++++------- Testing/ImageTests/ImageTests.cpp | 65 +++++++++++++++++++++++++++ 2 files changed, 123 insertions(+), 16 deletions(-) diff --git a/Libs/Image/ImageRegistration.cpp b/Libs/Image/ImageRegistration.cpp index 5a058110604..8100583edba 100644 --- a/Libs/Image/ImageRegistration.cpp +++ b/Libs/Image/ImageRegistration.cpp @@ -68,6 +68,11 @@ constexpr double MAXIMUM_STEP_FRACTION = 0.1; // as having lost the shape it was registering. See OverlapGuardCommand. constexpr double MINIMUM_OVERLAP_FRACTION = 0.5; +// How many points the overlap guard maps each iteration to see how much of the fixed image is still +// landing on the moving one. Enough to measure a fraction to a percent or so, few enough that the +// mapping costs nothing beside the metric evaluation it is watching over. +constexpr unsigned int OVERLAP_SAMPLE_TARGET = 512; + // The range of size differences the affine stage will be seeded with. Shapes of the same anatomy do // not differ by more than this, so an estimate outside it is a sign that the images are not what the // estimate assumes -- and the seed is only ever a starting point worth having, never one worth @@ -175,10 +180,18 @@ ScaleSeed estimate_scale_seed(const ImageType* fixed, const ImageType* moving) { /// that wanders into it stays, and the registration comes back mapping every point far outside the /// image it was supposed to land on. /// -/// How many samples the metric still had to work with is the tell. When that collapses, stop: the -/// optimizer holds on to the best parameters it found, which are the ones from before the escape, -/// and the remaining levels carry on from there. Stopping a level that legitimately needed to give -/// up half its overlap only costs some convergence, which is much the cheaper mistake. +/// How much of the fixed image is still landing on the moving one is the tell. This maps its own +/// fixed set of points through the transform to measure that, rather than reading the metric's count +/// of the samples it last found usable: with a line search that count belongs to whichever trial step +/// the search happened to evaluate last, not to the step actually taken, so it can condemn a step +/// that was fine and, worse, pass a step that was not. A few hundred points cost nothing against the +/// metric evaluation this is watching over, and being the same points every time makes the fractions +/// comparable across iterations and across levels of the pyramid. +/// +/// When the overlap collapses, stop: the optimizer holds on to the best parameters it found, which +/// are the ones from before the escape, and the remaining levels carry on from there. Stopping a +/// level that legitimately needed to give up half its overlap only costs some convergence, which is +/// much the cheaper mistake. template class OverlapGuardCommand : public itk::Command { public: @@ -186,7 +199,24 @@ class OverlapGuardCommand : public itk::Command { using Pointer = itk::SmartPointer; itkNewMacro(Self); - void set_metric(const TMetric* metric) { metric_ = metric; } + //! the transform to watch, and the two images whose overlap is being measured + void set_images(const TMetric* metric, const ImageType* fixed, const ImageType* moving) { + metric_ = metric; + moving_ = moving; + + const auto region = fixed->GetBufferedRegion(); + const auto stride = + std::max(itk::SizeValueType{1}, itk::SizeValueType{region.GetNumberOfPixels() / OVERLAP_SAMPLE_TARGET}); + itk::ImageRegionConstIteratorWithIndex it(fixed, region); + itk::SizeValueType counter = 0; + for (it.GoToBegin(); !it.IsAtEnd(); ++it, ++counter) { + if (counter % stride == 0) { + ImageType::PointType point; + fixed->TransformIndexToPhysicalPoint(it.GetIndex(), point); + samples_.push_back(point); + } + } + } void Execute(itk::Object* caller, const itk::EventObject& event) override { Execute(const_cast(caller), event); @@ -194,7 +224,7 @@ class OverlapGuardCommand : public itk::Command { void Execute(const itk::Object* caller, const itk::EventObject& event) override { auto* optimizer = const_cast(dynamic_cast(caller)); - if (!optimizer || !metric_) { + if (!optimizer || !metric_ || samples_.empty()) { return; } @@ -202,12 +232,21 @@ class OverlapGuardCommand : public itk::Command { return; } - // the count rises as the registration moves up the pyramid, so compare against the most this - // registration has managed rather than against a fixed number - const auto valid = metric_->GetNumberOfValidPoints(); - most_valid_ = std::max(most_valid_, valid); + const auto* transform = metric_->GetMovingTransform(); + if (!transform) { + return; + } + + itk::SizeValueType inside = 0; + for (const auto& sample : samples_) { + ImageType::IndexType index; + if (moving_->TransformPhysicalPointToIndex(transform->TransformPoint(sample), index)) { + inside++; + } + } - if (valid < static_cast(most_valid_ * MINIMUM_OVERLAP_FRACTION)) { + most_inside_ = std::max(most_inside_, inside); + if (inside < static_cast(most_inside_ * MINIMUM_OVERLAP_FRACTION)) { optimizer->StopOptimization(); } } @@ -217,14 +256,17 @@ class OverlapGuardCommand : public itk::Command { private: const TMetric* metric_{nullptr}; - itk::SizeValueType most_valid_{0}; + const ImageType* moving_{nullptr}; + std::vector samples_; + itk::SizeValueType most_inside_{0}; }; //--------------------------------------------------------------------------- /// Build an optimizer with the scale estimator wired to the given metric, so that rotation, /// translation and scaling parameters all step by a comparable physical distance. template -OptimizerType::Pointer make_optimizer(TMetric* metric, unsigned int iterations, const ImageType* fixed) { +OptimizerType::Pointer make_optimizer(TMetric* metric, unsigned int iterations, const ImageType* fixed, + const ImageType* moving) { using ScalesEstimatorType = itk::RegistrationParameterScalesFromPhysicalShift; auto scales_estimator = ScalesEstimatorType::New(); scales_estimator->SetMetric(metric); @@ -253,7 +295,7 @@ OptimizerType::Pointer make_optimizer(TMetric* metric, unsigned int iterations, optimizer->SetReturnBestParametersAndValue(true); auto guard = OverlapGuardCommand::New(); - guard->set_metric(metric); + guard->set_images(metric, fixed, moving); optimizer->AddObserver(itk::IterationEvent(), guard); return optimizer; @@ -389,7 +431,7 @@ void ImageRegistration::Impl::run_rigid() { registration->SetMetric(metric); registration->SetInitialTransform(rigid); registration->InPlaceOn(); - registration->SetOptimizer(make_optimizer(metric.GetPointer(), linear_iterations.front(), fixed)); + registration->SetOptimizer(make_optimizer(metric.GetPointer(), linear_iterations.front(), fixed, moving)); apply_multi_resolution_schedule(registration.GetPointer(), effective_linear_shrink_factors, linear_smoothing_sigmas); auto level_command = LevelIterationCommand::New(); @@ -425,7 +467,7 @@ void ImageRegistration::Impl::run_affine() { registration->SetMovingInitialTransform(composite); registration->SetInitialTransform(affine); registration->InPlaceOn(); - registration->SetOptimizer(make_optimizer(metric.GetPointer(), linear_iterations.front(), fixed)); + registration->SetOptimizer(make_optimizer(metric.GetPointer(), linear_iterations.front(), fixed, moving)); apply_multi_resolution_schedule(registration.GetPointer(), effective_linear_shrink_factors, linear_smoothing_sigmas); auto level_command = LevelIterationCommand::New(); diff --git a/Testing/ImageTests/ImageTests.cpp b/Testing/ImageTests/ImageTests.cpp index beb9630b693..b39b5d49b1d 100644 --- a/Testing/ImageTests/ImageTests.cpp +++ b/Testing/ImageTests/ImageTests.cpp @@ -1232,6 +1232,71 @@ TEST(ImageTests, registrationTranslationTest) { } } +TEST(ImageTests, registrationScaleTest) { + Image fixed_dt(std::string(TEST_DATA_DIR) + "/ellipsoid_00.DT.nrrd"); + // Big enough that the affine stage does not find it unaided: left to search from a pose the rigid + // stage had no freedom to scale, it returns a scale of 1 and maps the shape onto part of its + // target. Smaller differences it manages by itself and would not test anything. + const double scale = 2.0; + + // the same shape at a different size: stretching the grid and the values alike is what keeps the + // result a valid distance transform rather than one whose values no longer match its geometry + Image moving_dt(fixed_dt); + const auto spacing = fixed_dt.spacing(); + moving_dt.setSpacing(Vector3({spacing[0] * scale, spacing[1] * scale, spacing[2] * scale})); + moving_dt *= scale; + + // the linear stages are what has to find a size difference, so leave the deformable one out of it + ImageRegistration registration; + registration.set_transform_type(ImageRegistration::TransformType::Affine); + registration.run(ImageRegistration::make_registration_image(fixed_dt), + ImageRegistration::make_registration_image(moving_dt)); + + // How far apart mapped points end up is what says whether the size difference was found. The + // bound is loose on purpose: the point is that most of the difference is recovered rather than + // none of it, and pinning the figure down would only make this brittle. + const std::vector points = {Point3({-10.0, 0.0, 0.0}), Point3({10.0, 0.0, 0.0}), + Point3({0.0, 0.0, -12.0}), Point3({0.0, 0.0, 12.0})}; + auto transformed = registration.transform_points(points); + + for (size_t i = 0; i < points.size(); i += 2) { + const double before = points[i].EuclideanDistanceTo(points[i + 1]); + const double after = transformed[i].EuclideanDistanceTo(transformed[i + 1]); + ASSERT_GT(after / before, 1.0 + (scale - 1.0) / 2.0) + << "the affine stage recovered little of a " << scale << "x size difference"; + } +} + +//--------------------------------------------------------------------------- +TEST(ImageTests, registrationSettingsDescriptionTest) { + // Anything that changes the transform has to change the description, because that is what tells a + // caller holding a saved transform whether it is still the transform it would get today. A + // setting missing from it is a stale transform served as a fresh one. + ImageRegistration reference; + const std::string description = reference.settings_description(); + + auto changed = [&](const std::function& change) { + ImageRegistration registration; + change(registration); + return registration.settings_description() != description; + }; + + ASSERT_FALSE(changed([](ImageRegistration&) {})) << "identical settings must describe identically"; + + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_transform_type(ImageRegistration::TransformType::Rigid); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_gradient_step(0.5); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_update_field_variance(1.0); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_total_field_variance(1.0); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_iterations({1, 2, 3}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_shrink_factors({8, 4, 2}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_smoothing_sigmas({4.0, 2.0, 1.0}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_linear_iterations({1, 2, 3, 4}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_linear_shrink_factors({8, 6, 4, 2}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_linear_smoothing_sigmas({4.0, 3.0, 2.0, 1.0}); })); + ASSERT_TRUE(changed([](ImageRegistration& r) { r.set_correlation_radius(2); })); +} + +//--------------------------------------------------------------------------- TEST(ImageTests, registrationParticleTransferTest) { Image reference_dt(std::string(TEST_DATA_DIR) + "/ellipsoid_00.DT.nrrd"); Image target_dt(std::string(TEST_DATA_DIR) + "/ellipsoid_01.DT.nrrd");