I have written the following code, and obtained the correct maximum value from the FFT Normalized Cross Correlation tools in SimpleITK.
I am now wondering, how I can get this maximum value location(x,y) within the image?
Please find here a functioning version of FFT based translation initialization, based on the code I had previously shared:
It uses the uses a sequence of the RegionalMaxima, ConnectedComponents and LabelStatistics image filter to locate the maximum in the correlation image, and then produces a translation matrix.
I was able to successfully get this to work, but can a few questions about the FFT output/results.
Are the calculated peak locations w.r.t. the FFT output image or the fixed image?
(ie. if the peak location is 25,25 in space, would that be in regard to the output image or the fixed?)
The moving image is a small section of an undeformed image. I would like to calculate the difference between its original location in the undeformed image and its new “matched” location in the deformed/fixed image.
However, I am not sure how to do this in sitk because I am still confused about the different spatial domains that are used.
Will I need to somehow resample/resize the smaller undeformed full image into the deformed image space, and find the original subimage center point in the newly resized image, then calculate the difference between that point and the FFT peak point results?
In the code I liked to above, a translation transform is returned which maps physical points from the fixed image to the moving image. The details of the calculation are performed in the code.
If you have any specific questions about that code please let me know.
If you are confused about the spacial domains, I suggest reading:
So, for example, in a 2D case, if the translation calculated by FFT is [-50,80], the translation transform being returned is just applying this singular translation pair to any point (as if it is linearly changing within the image)?
I have been using some similar code as you attached but for 3D Image data.
The code you referenced finds the peak point using the bounding box location and converting it to physical space. How would this be extendable to 3 dimensions?
This is what I have right now, however, the resulting peak locations are being calculated as negative values which would not make sense in this case.
Thank you in advance!
Code:
sitk::Image out = sitk::FFTNormalizedCorrelation(fft_fixed, fft_moving);
out = sitk::SmoothingRecursiveGaussian(out);
//Find maximum
sitk::MinimumMaximumImageFilter minMax;
minMax.Execute(out);
double max = minMax.GetMaximum();
//Index max location
sitk::Image cc = sitk::ConnectedComponent(sitk::RegionalMaxima(out, 0.0, 1.0, true));
sitk::LabelStatisticsImageFilter stats;
stats.Execute(out, cc);
auto labels = stats.GetLabels();
std::vector<std::pair<double, int64_t>> pairs;
for (int i = 0; i < labels.size(); i++)
{
pairs.push_back(std::make_pair(stats.GetMean(labels[i]), labels[i]));
}
std::sort(pairs.begin(), pairs.end());
auto bb = stats.GetBoundingBox(pairs[pairs.size() - 1].second);
I would suggest you break this code down into separate testable steps. At each step the output can be visually inspected, and the quatitative results verified.
For example can procedurally create a set of label and intensity image to test different conditions. This will verify that the computation is correct and as you expect.
