Issue with Point Set Metrics in antsRegistration

[jjbouza]

Thank you all for the great work on this package! I am trying to use landmark points to guide the registration of some abdominal MR volumes and have encountered an unexpected issue. Here is the command I am trying to run:

antsRegistration -d 3 --metric ICP[fixed.txt,moving.txt,1] --metric MI[fixed.nii,moving.nii,1] --transform SyN[0.1,0,0] --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence [1000x100x100,1e-4,4] -v

The two volumes are already affine registered. I've attached my fixed.txt and moving.txt files (you should be able to reproduce with any fixed.nii and moving.nii volumes).

The issue with the command above is that the metricValue never decreases and the result of the registration is an identity transformation. When the --transform is changed to any other deformable transformation method (e.g. Exponential) I get a segmentation fault.

ANTs version 2.3.5, compiled from source on Ubuntu 20.04. I've reproduced the issue on a separate machine as well.

I appreciate the help.
moving.txt
fixed.txt

[cookpa]

The points are defined in the physical space of the images, so they are not independent of the fixed and moving images.

The points must be defined in ITK coordinates (LPS, like DICOM). You can load your images into ITK-SNAP and then see whether your points are sensible (Tools -> Layer Inspector, select Info tab).

[jjbouza]

Thanks @cookpa, I was using RAS coordinates. I mapped my coordinates to LPS and checked with itk-snap that they were sensible for my images. I am still seeing the same issue though, with the metricValue staying constant.

[jjbouza]

Here is a minimal example of the issue: https://drive.google.com/drive/folders/1aMS0gJFQsDZU9AKWKaTnrpxaLhqffa6z?usp=sharing

Running antsRegistration -d 3 -o "[warped,warped.nii]" --metric "ICP[fixed.txt,moving.txt,1.0]" --metric "MI[fixed.nii,moving.nii.gz,1.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

Gives:

XXDIAGNOSTIC,Iteration,metricValue,convergenceValue,ITERATION_TIME_INDEX,SINCE_LAST
 1DIAGNOSTIC,     1, 4.029707796358e+00, 1.797693134862e+308, 1.5200e-01, 1.5200e-01,
 1DIAGNOSTIC,     2, 4.029707796358e+00, 1.797693134862e+308, 2.0391e-01, 5.1909e-02,
 1DIAGNOSTIC,     3, 4.029707796358e+00, 1.797693134862e+308, 2.5408e-01, 5.0177e-02,
 1DIAGNOSTIC,     4, 4.029707796358e+00, -1.700143779632e-05, 3.0570e-01, 5.1615e-02,
XXDIAGNOSTIC,Iteration,metricValue,convergenceValue,ITERATION_TIME_INDEX,SINCE_LAST
 1DIAGNOSTIC,     1, 4.029707796358e+00, 1.797693134862e+308, 8.2418e-01, 5.1848e-01,
 1DIAGNOSTIC,     2, 4.029707796358e+00, 1.797693134862e+308, 1.1988e+00, 3.7461e-01,
 1DIAGNOSTIC,     3, 4.029707796358e+00, 1.797693134862e+308, 1.5587e+00, 3.5994e-01,
 1DIAGNOSTIC,     4, 4.029707796358e+00, -1.700143779632e-05, 1.9442e+00, 3.8551e-01,
XXDIAGNOSTIC,Iteration,metricValue,convergenceValue,ITERATION_TIME_INDEX,SINCE_LAST
 1DIAGNOSTIC,     1, 4.029707796358e+00, 1.797693134862e+308, 6.2566e+00, 4.3123e+00,
 1DIAGNOSTIC,     2, 4.029707796358e+00, 1.797693134862e+308, 1.0251e+01, 3.9944e+00,
 1DIAGNOSTIC,     3, 4.029707796358e+00, 1.797693134862e+308, 1.3757e+01, 3.5063e+00,
 1DIAGNOSTIC,     4, 4.029707796358e+00, -1.700143779632e-05, 1.7005e+01, 3.2479e+00,
  Elapsed time (stage 0): 17.2794

[ntustison]

Have you looked at this example repository?

[jjbouza]

@ntustison Yes, I did. Correct me if I'm wrong, but I don't think that repository has an example of registering point sets specified as space delimited .txt files. E.g. antsRegistrationCommandDiffeomorphismTest.sh in that repo has fixedPoints='data/pointsLowerLeft.nii.gz' which I assume is a NIFTI image where landmark points are represented as disks with some non-zero radius. I would rather pass the landmark points in a txt file if possible.

[ntustison]

You're right but you can convert those labeled images to .csv point set files using ImageMath ... LabelStats ..., ensure that it works as expected, and then map it to your own problem.

[cookpa]

@ntustison is there a minimum number of points for the metric to work effectively? The files above only have two points

[ntustison]

I hadn't noticed that. Thanks for pointing that out. Two points should work, I think, but I wouldn't use it as a case for debugging.

[jjbouza]

@ntustison is there a minimum number of points for the metric to work effectively? The files above only have two points

This is a minimal example, I tried with up to 6 point correspondences with the same issue.

[jjbouza]

You're right but you can convert those labeled images to .csv point set files using ImageMath ... LabelStats ..., ensure that it works as expected, and then map it to your own problem.

I confirmed that the following commands reproduce the issue with the "points" data in the chicken repository:

git clone https://github.com/stnava/chicken.git
cd chicken

# Convert *.nii.gz files to *.csv
ImageMath 2 data/pointsLowerLeft.csv LabelStats data/pointsLowerLeft.nii.gz data/pointsLowerLeft.nii.gz
ImageMath 2 data/pointsUpperRight.csv LabelStats data/pointsUpperRight.nii.gz data/pointsUpperRight.nii.gz

# replace commas with spaces
cat data/pointsLowerLeft.csv | tr ',' ' ' > data/pointsLowerLeft.txt
cat data/pointsUpperRight.csv | tr ',' ' ' > data/pointsUpperRight.txt

# remove header
sed -i '1d' data/pointsLowerLeft.txt
sed -i '1d' data/pointsUpperRight.txt

# select cols
cut -d ' ' -f 1,2,3,5 data/pointsLowerLeft.txt > data/pointsLowerLeftFinal.txt
cut -d ' ' -f 1,2,3,5 data/pointsUpperRight.txt > data/pointsUpperRightFinal.txt

# run registration
antsRegistration -d 2 -o "[warped,warped.nii]" --metric "ICP[data/pointsLowerLeftFinal.txt,data/pointsUpperRightFinal.txt,1.0]" --metric "MI[data/pointsLowerLeft.nii.gz,data/pointsUpperRight.nii.gz,1.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

[ntustison]

In my suggestion above, my assumption was that you'd do a bit more digging on your end. The original chicken example works so the key question is "at what stage in making it more like your problem does it not work?" For example, is it the ICP metric? Is it the .txt conversion? We're busy with our own projects, etc. and so we rely on users to do as much as they can before we take on doing the investigation ourselves.

[jjbouza]

Thanks @ntustison, I can definitely offer some more details (and I appreciate the continued help, I know this is a side-project). I have dug into this for a couple of days but I'm trying not to present too much information at once. Here are some things I've noticed:

• The issue is not specific to ICP, it happens with all point-set metrics.
• I think the format of the *.txt file is correct (I'm basing it of a few previous Github issues), but I would appreciate if you could review it.
• When the order of the MI and ICP metrics is switched, the registration metric value actually does decrease, but I am skeptical that the ICP metric is being used in the registration, since the landmark alignment after registration is just as poor as only using the MI metric. This point is also confusing to me because I would not expect that switching the order of the MI and ICP metrics would change the registration results (beyond the expected small random deviations in registration performance).
• When the ICP metric is used as the only metric, the registration does not run at all (segfault).

My conclusion from this is that the *.txt landmark point metric losses are not being used to update the transformation parameters internally. It is unclear to me why this would be the case. The original chicken example works because it uses *.nii.gz files with "landmark disks", but I have not been able to find a working example with *.txt landmark points.

If we can figure this out together I would be happy to contribute some documentation and examples for running registration with *.txt landmark points.

[jjbouza]

Also relevant: The same problem was touched on in this thread but never resolved.

[ntustison]

I don't know what "landmark disks" are but the input data structure of the point set metric is the same for images as it is for .txt files.

[jjbouza]

By "landmark disks" I just mean that e.g. data/pointsLowerLeft.nii.gz in the chicken repository represents a "landmark" as a disk shape region of pixels with non-zero radius. I agree that internally they are both mapped to the same data structure.

Maybe a good place to start debugging this is a simpler example: what is the expected behavior of registration with only ICP metric?

antsRegistration -d 2 -o "warped" --metric "ICP[data/pointsLowerLeftFinal.txt,data,pointsUpperRightFinal.txt,1.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

This results in a segfault for me:

All_Command_lines_OK
Using double precision for computations.
  number of levels = 3
  fixed labeled point set: data/pointsLowerLeftFinal.txt
  moving labeled point set: data/pointsUpperRightFinal.txt
Dimension = 2
Number of stages = 1
Use Histogram Matching false
Winsorize image intensities false
Lower quantile = 0
Upper quantile = 1
Stage 1 State
   Point Set Metric = ICP
     Fixed labeled point set = PointSet (0x55a307d760d0)
  RTTI typeinfo:   itk::PointSet<unsigned int, 2u, itk::DefaultStaticMeshTraits<unsigned int, 2u, 2u, float, float, unsigned int> >
  Reference Count: 2
  Modified Time: 355
  Debug: Off
  Object Name:
  Observers:
    none
  Source: (none)
  Source output name: (none)
  Release Data: Off
  Data Released: False
  Global Release Data: Off
  PipelineMTime: 0
  UpdateMTime: 316
  RealTimeStamp: 0 seconds
  Number Of Points: 4
  Requested Number Of Regions: 1
  Requested Region: 0
  Buffered Region: -1
  Maximum Number Of Regions: 1
  Point Data Container pointer: 0x55a307bf3bf0
  Size of Point Data Container: 4

     Moving labeled point set = PointSet (0x55a307d77a30)
  RTTI typeinfo:   itk::PointSet<unsigned int, 2u, itk::DefaultStaticMeshTraits<unsigned int, 2u, 2u, float, float, unsigned int> >
  Reference Count: 2
  Modified Time: 356
  Debug: Off
  Object Name:
  Observers:
    none
  Source: (none)
  Source output name: (none)
  Release Data: Off
  Data Released: False
  Global Release Data: Off
  PipelineMTime: 0
  UpdateMTime: 353
  RealTimeStamp: 0 seconds
  Number Of Points: 4
  Requested Number Of Regions: 1
  Requested Region: 0
  Buffered Region: -1
  Maximum Number Of Regions: 1
  Point Data Container pointer: 0x55a307d720a0
  Size of Point Data Container: 4

     Weighting = 1
     Use only boundary points = false
     Point set sigma = 1
     Evaluation K neighborhood = 50
     Alpha = 1.1
     Use anisotropic covariances = false
     Sampling percentage = 1
   Transform = SyN
     Gradient step = 0.1
     Update field sigma (voxel space) = 0
     Total field sigma (voxel space) = 0
     Update field time sigma = 0
     Total field time sigma  = 0
     Number of time indices = 0
     Number of time point samples = 0
Registration using 1 total stages.

Stage 0
  iterations = 1000x100x100
  convergence threshold = 0.0001
  convergence window size = 4
  number of levels = 3
  using the ICP metric (weight = 1)
[1]    4420 segmentation fault  antsRegistration -d 2 -o "[warped,warped.nii]" --metric  --transform   4x2x1

Does ANTs need the header info in NIFTI files to run registration?

[cookpa]

I think you need to have an image metric. It's needed to define physical space of the inputs and outputs. I'm not familiar with all the point set data types, but a text file is definitely not enough to initialize the registration. In any case, I don't see how you could solve for a dense warp field from a handful of distant points.

It might be more theoretically reasonable to define an affine transform with points only, but I don't know if antsRegistration can do that.

You could try using Mattes with a small weight to approximate a solution that relies more on the point set metric.

[jjbouza]

Thanks @cookpa. So I think we can safely assume --metric ICP[fixed.txt,moving.txt,1.0] does not work because there is insufficient information to define the physical space.

So then we have two other (possibly related) issues:

• --metric ICP[fixed.txt,moving.txt,1] --metric MI[fixed.nii,moving.nii,1] does not segfault, but results in a constant metricValue: Issue with Point Set Metrics in antsRegistration #1335 (comment). Maybe antsRegistration will still try to pull the physical space parameters from ICP since it is the first metric, but will fail more gracefully?
• --metric MI[moving.nii,fixed.nii,1] --metric ICP[fixed.txt,moving.txt,1] works, but it appears that the ICP metric has no effect on the generated transformation. I infer this from the fact that assigning a large weight to ICP --metric MI[moving.nii,fixed.nii,1] --metric ICP[fixed.txt,moving.txt,10000] results in the same transformation. Also, completely removing the ICP metric --metric MI[moving.nii,fixed.nii,1] results in the same transformation again. In both cases, the generated transformations are not correct.

[ntustison]

I believe I wrote it such that one could specify a mask to provide the necessary physical domain without having to use an image-specific metric. Also, you might want to explore this issue using MeasureImageSimilarity to isolate any metric issues from the rest of the registration machinery.

[cookpa]

Good to know about the mask. I'll look into this when I get time and see about updating the tutorials and documentation

[jjbouza]

The mask comment helps a lot. The following command fixes the issue with the segfault from here:
antsRegistration -d 2 -o "warped" --masks "[data/pointsLowerLeft.nii.gz,data/pointsUpperRight.nii.gz]" --metric "ICP[data/pointsLowerLeftFinal.txt,data,pointsUpperRightFinal.txt,1.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

Now the registration appears to be working for the ICP only case. I'll experiment with adding an image similarity metric as well and get back to you.

[jjbouza]

I generated a binary mask image in the same physical space as data/pointsLowerLeft.nii.gz and data/pointsUpperRight.nii.gz with all pixels set to value 1. I confirmed that ICP only registration still worked with this mask:
antsRegistration -d 2 -o "warped" --masks "[data/mask.nii.gz,data/mask.nii.gz]" --metric "ICP[data/pointsLowerLeftFinal.txt,data,pointsUpperRightFinal.txt,1.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

Next I tried adding --MI metric with a weight of zero:
antsRegistration -d 2 -o "warped" --masks "[data/mask.nii.gz,data/mask.nii.gz]" --metric "ICP[data/pointsLowerLeftFinal.txt,data,pointsUpperRightFinal.txt,1.0]" --metric "MI[data/pointsLowerLeft.nii.gz,data/pointsUpperRight.nii.gz,0.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

I expected to get a similar result to the previous command, but instead the metricValue is again constant throughout registration. Interestingly, the generated warp field is not an identity deformation, but it is incorrect.

The same issue occurs when using *.nii.gz inputs to the ICP metric, i.e. the following command gives the same incorrect result:
antsRegistration -d 2 -o "warped" --masks "[data/mask.nii.gz,data/mask.nii.gz]" --metric "ICP[data/pointsLowerLeft.nii.gz,data,pointsUpperRight.nii.gz,1.0]" --metric "MI[data/pointsLowerLeft.nii.gz,data/pointsUpperRight.nii.gz,0.0]" --transform "SyN[0.1,0,0]" --shrink-factors 4x2x1 --smoothing-sigmas 3x1x0 --convergence "[1000x100x100,1e-4,4]" -v

So the problem is independent of the input types passed to --metric ICP. It seems there is an internal issue when using both an image similarity metric and a point set metric. Any ideas for investigating further?

[cookpa]

I see what you're saying, the metric values don't seem to change. I will look more into this when I get time.

Is there a reason you're using --transform "SyN[0.1,0,0]"? I think not smoothing the update field at all could cause convergence issues as well, though not as extreme as this. I'd stick with SyN[0.1,3,0] for testing.

[jjbouza]

Thanks @cookpa, I spent some time last week going through the source code but could not find where the issue was coming from. One more note that might be useful for debugging: the issue does not appear when using affine or rigid registration (I managed to get that to work), but does appear for all deformable transformation models.

--transform "SyN[0.1,0,0]" was used for testing purposes to eliminate any potential reasons that the landmarks would not be aligned after registration. I tried with update field smoothing greater than 0 as well.

[ntustison]

For testing purposes, I would use "BSplineSyN" as it's much more appropriate framework for both voxel and landmark contributions.

[cookpa]

I've not made any breakthroughs here but as we've been dealing with metric scales issues lately, I began to wonder if the scales estimation might also be relevant to this issue. I noticed scales are set by the first metric. It might not be the only time the first metric is important, but I thought it might be a lead to explain why the ordering of metrics matters in the experiments above

ANTs/Examples/itkantsRegistrationHelper.hxx

Lines 1315 to 1321 in ec91afc

	// There's a scale issue here. Currently we are using the first metric to estimate the
	// scales but we might need to change this.
	typedef itk::RegistrationParameterScalesFromPhysicalShift<ObjectMetricType> ScalesEstimatorType;
	typename ScalesEstimatorType::Pointer scalesEstimator = ScalesEstimatorType::New();
	scalesEstimator->SetMetric(singleMetric);
	scalesEstimator->SetTransformForward(true);