ImageRegistrationMethodv4 VirtualImage with SparseSampling

blowekamp · July 27, 2018, 1:05pm

When using the ImageRegistrationMethod with a virtual image and a sampling strategy, there is a serious problem with the spaces.

->The ImageRegistrationMethod generate a sampled point set by iterating on the virtual image:

github.com

InsightSoftwareConsortium/ITK/blob/master/Modules/Registration/RegistrationMethodsv4/include/itkImageRegistrationMethodv4.hxx#L954-L1035


      
          if (m_ReseedIterator)
          {
            randomizer->SetSeed();
          }
          else
          {
            randomizer->SetSeed(m_CurrentRandomSeed++);
          }
          
          
          unsigned long index = 0;
          
          switch (this->m_MetricSamplingStrategy)
          {
            case MetricSamplingStrategyEnum::REGULAR:
            {
              const auto sampleCount =
                static_cast<unsigned long>(std::ceil(1.0 / this->m_MetricSamplingPercentagePerLevel[this->m_CurrentLevel]));
              unsigned long count =
                sampleCount; // Start at sampleCount to keep behavior backwards identical, using first element.

This file has been truncated. show original

This is then set to the metric by ImageToImageMetricv4::SetFixedSampledPointSet. The Metric expects the samples to be in the fixed image domain!

→ The metric then transforms these points from the fixed image domain to the virtual image domain :

github.com

InsightSoftwareConsortium/ITK/blob/master/Modules/Registration/Metricsv4/include/itkImageToImageMetricv4.hxx#L520-L567


      
          
          template <typename TFixedImage,
                    typename TMovingImage,
                    typename TVirtualImage,
                    typename TInternalComputationValueType,
                    typename TMetricTraits>
          void
          ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>::
            SetMaximumNumberOfWorkUnits(const ThreadIdType number)
          {
            if (number != this->m_SparseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads())
            {
              this->m_SparseGetValueAndDerivativeThreader->SetMaximumNumberOfThreads(number);
              this->m_SparseGetValueAndDerivativeThreader->SetNumberOfWorkUnits(number);
              this->Modified();
            }
            if (number != this->m_DenseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads())
            {
              this->m_DenseGetValueAndDerivativeThreader->SetMaximumNumberOfThreads(number);
              this->m_DenseGetValueAndDerivativeThreader->SetNumberOfWorkUnits(number);

This file has been truncated. show original

I have a patch which has a hack for a solution for how I was using it:
http://review.source.kitware.com/#/c/23593/

Here is the commit message:

For using a sampled point set for a metric in the
ImageRegistrationMethodv4, the point set was generated by iterating
over the virtual domain. This matches the way the dense sampling is
done. However, this is set to the metrics as the fixed point set,
which the metric then transform from the fixed image domain to the
virtual image domain, essentially applying the initial fixed image
transform twice.

The SetFixedSampledPointSet is widely used. and the
“UseFixedSampledPointSet” ivar is checked to determine if the
“VirtualSampledPointSet” can be used. The FixedSampledPointSet is only
used by the initialization of the metric, as the evaluation of the
metrics access the virtual sampled point set.

It could be made clear by a “UseSampledPointSet” variable instead of
the “UsedFixedSampledPointSet”, but this would be an API change.

When the VirtualSampledPointSet is set should it be transformed to the
FixedSampledPointSet? This data is not used, but would help keep the
state of the class consistent.

Steps to move forward:

Agreement this is a bug
Determine compatibility requirements.
->Need a way to clearly set fixed or virtual point set, but currently UsedFixedSampledPointSet is use to determine if the transformed virtual point set is used, while the original point set does not appear to be used anywhere. Keeping UseFixedSamplesPointSet will be very confusing, misleading, and a bad name for the states of the variables!

blowekamp · July 27, 2018, 1:08pm

@hjmjohnson

hjmjohnson · July 27, 2018, 1:28pm

Can you reach out to Nick Tustison first?

ntustison · August 5, 2018, 4:22pm

Hey Brad,

Sorry for not getting back to you earlier.

I looked at what you did and from what I can tell, your suggested
changes make sense. It sounds like it was an oversight on our
part in not maintaining a strict separation between the fixed and
virtual domains. That might be explained as I wrote the image
registration methods while Michael wrote the image to image
metrics. Brian watched over both of our contributions but it’s easy
to see how something like this could slip through.

When you outline the problem, I’m guessing it hasn’t translated into
incorrect registration solutions because the transformation between
the fixed and virtual domains is, by default, the identity transform.

Anyway, you mention your fix having a “hack”. Is there something
more I can do to help you with this solution?

Nick

blowekamp · August 6, 2018, 1:58pm

Nick,

Good to head from you!

I have not looked too closely at the Syn methods. I thought they used the virtual domain more especially when used with group wise registration. The deformation fields use dense sampling, so no issue there, but the BSpline transforms may use a sample point sets which would be problematic. Is this virtual image commonly used any where in ANTs?

Regarding correctness when this issue is encountered, because the are just the sample points the resulting transform is still correct. The problem would be that the fixed image transform is applied twice, so if it small it might not be noticed but the expected sampling domain would not be as described. When the fixed transform is large all of the point will be mapped out of the fixed image, and exceptions are generated.

In my “steps to move forward”, it looks like there is agreement there is a bug, (#1), now we need to figure out (#2) what IVARs or setters are needed.

I referred to my patch as a hack be cause it was quickly done with out much testing. As mentioned above anding a state variable to deterring the if VirtualSampledPointSet or FixedSampledPointSet IVAR is the one set by the user is needed but with the current UseFixedSampledPointSet this will be confusing. So to have sensible setter name I’ll proposer the following

Rename:
UseFixedSampledPointSet -> UseSampledPointSet
- UseFixedSampledPointSet would be marked as deprecated and modified to set/get the new m_SampledPointSet IVAR. Unfortunately m_UseFixedSampledPointSet is protected and used by derived classes, so there is there is some “internal/derived” compatibility issues because the class has protected access for it’s IVARs.
Add:
UseVirtualSampledPointSet

This will default to false, so that the FixedSampledPointSet is used by default for compatibility.

ntustison · August 6, 2018, 4:36pm

This sounds good to me.

I should point out that all the virtual domain stuff is currently a placeholder for potential future work. The virtual domain and fixed image domain are basically synonymous in the current ITKv4 framework such that the transform between them is always identity. This includes all the special image registration methods such as the SyN and other diffeomorphic variants. Although we tried to keep IVAR nomenclature consistent with this potential future work, mistakes could’ve easily been made.

Nick

blowekamp · August 6, 2018, 6:08pm

This patch is also more efficient since it remove a copying of the point sent, so there really is no down side to it.

I encountered this problem while working with microscopy images that have a large number of pixels, but small pixel space. The question came about what would be the best pixel spacing to perform image registration with respects to numerics. We can hack at the image’s spacing directly, or see if a transform to a standard/normalized virtual image domain works.