TransformPhysicalPointToIndex that remains inside the boundaries of the image

Hello everyone,

I am working with ITK images and I am using TransformPhysicalPointToIndex a lot!.
My physical space bounds are quite strict, so I wanna be sure that this function computes accurately indices of points that are inside an image or not.

I have a problem that is related to the fact that this function is using rounding.

I have no problem doing rounding for the interior points of the image.
But when it comes to the boundary, if a point is even slightly outside of the image I don’t want this function to return me an index that is part of the image, when geometrically speaking it’s not.

Is there a similar function like TransformPhysicalPointToBoundedIndex(point, index)
that implements the logic that I explained?

Hello @spyridon97,

Possibly use TransformPhysicalPointToContinuousIndex and then you can make the decision of how to deal with the returned index?

That could be a solution, but since I use this function in many places I would prefer something built-in that does the job by itself. if there is not something like that, do you think that ITK library would be open to such an addition?

I think we would be open. I needed such a thing around 10 years ago :smiley:

Ok then! I will try to implement it inside the ImageBase class and create a merge request for it.

I have created a pull request for it.

Hello,

Thank you for your contribution :clap:

I have left some comments on the PR regarding the ITK’s concept of indexing and bounds. That is for a 1-D image with a [0,10] region the valid continuous index range is [-0.5, 9.5]. So your method will need to correctly account for that.

Thanks for implementing your idea @spyridon97. But it turned out that what you wanted was different definition of inside, not just more strictly checked one. I am sorry your first patch did not have a happy ending :cry:

1 Like

It’s ok. Can you explain what you needed 10 years ago?

1 Like

I think that 10 years ago the pixel origin convention was not consistently applied within the ITK, and that was one of the goals of refactoring done in ITKv4. I was also prominently using VTK in my program too, and VTK had different definition of pixel origin. It took me a while to figure it all out, with occasional crashes due to some vertex going out of image by half a voxel or so.

I thought you were going to improve numerical precision of current inside/outside decision-making.

Although I can see the appeal of avoiding the problem by having a half a voxel wide safety margin, having it integrated into ITK would make it harder for people to grasp all the physical space and index space concepts. Not to mention that it is mostly being motivated by using a viewer with different definition of pixel origin.

Ok, now understand all the decision making behind that choice of the origin and the differences between vtk image and itk image.

So, the ITK library would like to have better precision when it comes to inside/outside decision-making.
and more precisely, this part of IsInside(const ContinuousIndex<TCoordRepType, VImageDimension> & index) const

      /* Note for NaN: test using negation of a positive test in order
       * to always evaluate to true (and thus return false) when index[i]
       * is NaN. The cast above to integer via RoundHalfIntegerUp will cast
       * NaN into a platform-dependent value (large negative, -1 or large
       * positive, empirically). Thus this test here is relied on
       * to 'catch' NaN's. */
      if (!(index[i] <= bound))
      {
        return false;
      }

I think that this could be done by utilizing exact-numeric-analysis principles. I will see what I can do regarding this issue.

Thanks a lot for explaining to me all these details.