Atropos stops early with "posterior probability" that doesn't depend on intensities data

stephen-q-le commented 1 month ago

Operating system and version

NixOS 24.05.3103.0c53b6b8c2a3 (Uakari)

CPU architecture

x86_64 (PC, Intel Mac, other Intel/AMD)

ANTs code version

ANTs Version: 2.5.1.dev1-gGIT-NOTFOUND Compiled: Jan 1 1980 00:00:00

ANTs installation type

Other (please specify below): NixOS package install

Summary of the problem

I have a brain mask such that the Atropos output looks like this, for various image intensities:

$ Atropos -v -d 3 -a <various> -x redacted-brainmask-map.nii.gz -i <various>[3] -c [5,0.0001] -m <various> -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]

Running Atropos for 3-dimensional images.

Progress: 
  Iteration 0 (of 5): posterior probability = 0 (annealing temperature = 1)
  Iteration 1 (of 5): posterior probability = 0.596973 (annealing temperature = 1)

Writing output:

This seems to be a bug, as neither of these make sense to me:

The "posterior probability" of Iteration 1 being the same for different intensity images, and initializations.
Atropos stopping before convergence. I encountered this issue when trying to segment a brain CT scan, and could see that the output was noticeably different when I forced more iterations with -c [5,0].

I don't understand why Atropos would ever stop if it hasn't reached the number of iterations limit and hasn't generated two "posterior probability" values close to each other, but I'm not sure that's a bug.

Commands to reproduce the problem.

Atropos -v -d 3 -a garbage_data-label.nii.gz -x redacted-brainmask-map.nii.gz -i random[3] -c [5,0.0001] -m [1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]
Atropos -v -d 3 -a garbage_data-label.nii.gz -x redacted-brainmask-map.nii.gz -i random[3] -c [5,0.0001] -m [0.1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]
# different intensities file, same behavior:
Atropos -v -d 3 -a garbage_data-vesselness.nii.gz -x redacted-brainmask-map.nii.gz -i random[3] -c [5,0.0001] -m [1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]
# k-means initialization, same behavior:
Atropos -v -d 3 -a garbage_data-vesselness.nii.gz -x redacted-brainmask-map.nii.gz -i kmeans[3] -c [5,0.0001] -m [1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]

# These commands are not included because they result in unrelated errors:
# Atropos -v -d 3 -a garbage_data-label.nii.gz -x redacted-brainmask-map.nii.gz -i kmeans[3] -c [5,0.0001] -m [1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz]
# Atropos -v -d 3 -a garbage_data-label.nii.gz -x redacted-brainmask-map.nii.gz -i kmeans[3] -c [5,0.0001] -m [1,1x1x1] -o [output/classified.nii.gz,output/posterior%02d.nii.gz]

Output of the command with verbose output.

full-log.txt

Data to reproduce the problem

redacted-brainmask-map.nii.gz garbage_data-label.nii.gz Note: I encountered this issue on a high-resolution CT scan, with a brain mask I created outside ANTs. To replicate the issue using data I'm willing to upload, I kept the same image dimensions and resolution, and replaced the intensities with hand-drawn garbage data. I made up another garbage data file "garbage_data-vesselness.nii.gz" but it's too big to upload.

ntustison commented 1 month ago

Not a bug. The posterior probability has converged and the final posterior is not dumped to the screen. It's a quirk of the command iterator.

stephen-q-le commented 1 month ago

You mean that the mean maximum posterior probability is 0.596973 both before and after Iteration 1, and that's how Atropos decides to stop? That makes sense, but still leaves the question of why the posterior probability doesn't depend on the input intensities image or on the initialization.

It also doesn't depend on the number of classes: Atropos -v -d 3 -a garbage_data-label.nii.gz -x redacted-brainmask-map.nii.gz -i random[6] -c [5,0.0001] -m [1,1x1x1] -k HistogramParzenWindows -o [output/classified.nii.gz,output/posterior%02d.nii.gz] 6classes.txt ... which also seems wrong to me: If the same image is segmented into 6 classes instead of 3 classes, and each of the 6 classes is actually used with roughly equal proportions, shouldn't the mean maximum posterior probability in the first iteration be lower? How can it be that this value is always 0.596973 for this particular resolution and mask shape?

ntustison commented 1 month ago

For starters, your example image is probably contributing to the problematic output. I would recommend doing some reading on this type of segmentation approach if you're really curious about how things are working "under the hood." I would start by reading the original Atropos paper and perhaps some of its references.

stephen-q-le commented 1 month ago

I computed mean maximum posterior (or at least, what I understand "We track convergence by summing up the maximum posterior probability at each site over the segmentation domain." and "mean maximum posterior probability over the region of interest" to mean) from the output images like this (compute_mean_maximum_posterior.py):

max_posterior = np.maximum.reduce(posteriors)
masked_max_posterior = max_posterior[mask_data > 0]
mean_max_posterior = np.mean(masked_max_posterior)

and the results are completely different from the 0.596973 number in the Atropos output, and different from each other: For the output of different Atropos runs (which all stop after reporting posterior probability = 0.596973) I get the numbers 0.9782532804268764, 0.9855073267560869, 0.9995354349536029.

I can't find the definition of GetCurrentPosteriorProbability() in this repository, can you tell me where that is defined?

When I run Atropos on a different brain CT I downloaded (CTBrain.nii.gz; I created a brain mask for it which is brainmask.nii.gz, it behaves as I expect: The posterior probability = values change when I switch the initialization from kmeans[3] to random[3], and Atropos takes more steps when the initialization is random[3]. In that case, the mean maximum posterior I compute is similar to the last posterior probability = reported by Atropos: 0.999974 vs 0.9999531303537921, and 0.975874 vs 0.9785478443309583.

Full script and output: script-compute-max-posterior.sh output-compute-max-posterior.txt

For starters, your example image is probably contributing to the problematic output

So yes, something about my input files is contributing to this, but I've now observed it on very different intensity images:

The brain CT I want to segment
A blurred version of that brain CT
garbage_data-label.nii.gz, some hand drawn squiggles, where each voxel intensity is a small integer
garbage_data-vesselness.nii.gz, a filter applied to those squiggles, such that the voxel intensities smoothly vary between 0 and 1 (but most are 0) I've also observed it with two brain masks; my original brain mask (where I always get posterior probability = 0.72874) and the modified brain mask redacted-brainmask-map.nii.gz.

cookpa commented 1 month ago

Get / Set methods are often defined via macro, so a text search for GetSomething will not work, but you can often find m_Something in the code.

In the case of m_CurrentPosteriorProbability, I think it's updated here:

https://github.com/ANTsX/ANTs/blob/9aebf5abddee8963e7073d6a1579e97da27ba811/ImageSegmentation/antsAtroposSegmentationImageFilter.hxx#L1148-L1150

ntustison commented 1 month ago

And note that although Atropos is a generic segmentation tool, in the brain, the underlying methodology is most often applied to T1-w MRI and not CT. In fact, depending on what segmentation output you're looking for in your CT image, Atropos might not be the best tool.

stephen-q-le commented 1 month ago

I think I figured out what's going wrong: On my system (both the NixOS installation and when I compile it myself), RealType is equal to float, not double. My mask has 28103827 voxels, and in AtroposSegmentationImageFilter::UpdateClassLabeling(), maxPosteriorSum is set to 1.67772e+07 because that's the the largest float I can get by repeatedly adding 1.

I could fix this loop but I don't know what else has similar issues; is there some easy way I can force RealType to be double? (And why wouldn't RealType be at least double in the first place?)

cookpa commented 1 month ago

Good find - I think this might not be just your OS, because I can reproduce the same thing on my Mac. It doesn't happen if I downsample your example images to 1mm resolution.

cookpa commented 1 month ago

Looks like something we would need to fix. We often use a consistent 32-bit RealType to conserve memory and reduce overhead from casting, which can be significant. But things that are summed over whole volumes should probably be double, what do you think @ntustison ?

ntustison commented 1 month ago

Sure, that shouldn't be a problem. You might even be able to simply swap out RealType in the header file without too many additional changes.

cookpa commented 1 month ago

Thanks for your persistence in tracking this down, @stephen-q-le . Let us know if there's further issues

ANTsX / ANTs