TransformIndexToPhysicalPoint. Shouldn't it take into account the starting index of the region?

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!

Origin point corresponds to index 0,0,0. I think this simplifies calculations with region of interest and pad filters, and makes it mentally easier for humans to reason about it. If there are any inconsistencies, they are bugs.

3 Likes

That’s it. I tend to forget that and think that Origin is the start pixel of the region. Thanks @dzenanz!

1 Like