This is the implementation of TransformToPhysicalPoint in itkImageBase:
/** Get a physical point (in the space which
* the origin and spacing information comes from)
* from a discrete index (in the index space)
*
* \sa Transform */
template <typename TCoordRep>
void
TransformIndexToPhysicalPoint(const IndexType & index, Point<TCoordRep, VImageDimension> & point) const
{
for (unsigned int i = 0; i < VImageDimension; ++i)
{
point[i] = this->m_Origin[i];
for (unsigned int j = 0; j < VImageDimension; ++j)
{
point[i] += m_IndexToPhysicalPoint[i][j] * index[j];
}
}
}
Shouldn’t that index be computed from the start index of the region?
An example:
// Allocate an image with a region with starting index {100,100,100} and size {4,4,4}
typename ITKImage::IndexType start;
typename ITKImage::SizeType size;
start.Fill(100);
size.Fill(4);
const typename ITKImage::RegionType region(start,size);
myITKImagePointer->SetRegions(region);
myITKImagePointer->Allocate();
myITKImagePointer->FillBuffer(0);
// And set a non-default Origin.
typename ITKImageType::PointType new_origin;
new_origin.Fill(-20);
myITKImagePointer->SetOrigin(new_origin);
typename ITKImage::IndexType start_index = myITKImagePointer->GetLargestPossibleRegion().GetIndex();
typename ITKImageType::PointType physical_point;
myITKImagePointer->TransformIndexToPhysicalPoint(start_index, physical_point);
Shouldn’t we expect physical_point to be equal to GetOrigin() because we used the first index of the region? This is happening right now only when start_index is equal to the default {0,0,0}.
In the implementation of the transform function is quite clear that the input index acts more like an offset than an index from a region. But I wonder if there is some other overloaded function that has this into account, or if I am missing something.
The code will behave as I would expect if:
template <typename TCoordRep>
void
TransformIndexToPhysicalPoint(const IndexType & index, Point<TCoordRep, VImageDimension> & point) const
{
IndexType normalized_index;
for (unsigned int i = 0; i < VImageDimension; ++i)
{
point[i] = this->m_Origin[i];
// Subtract the index of the largest possible region for regions with non-zero starting index.
normalized_index[k] = index[k] - this->GetLargestPossibleRegion().GetIndex()[k];
for (unsigned int j = 0; j < VImageDimension; ++j)
{
point[i] += m_IndexToPhysicalPoint[i][j] * normalized_index[j];
}
}
}
But this needs a region to be defined, and maybe more assumptions. Also not sure how many things this change will touch. I think only a few filters use a non-zero starting index (padFilters for example), and maybe that’s why this wasn’t spotted before.
What do you think? Thanks!