Hough Transform 2D Circles Image Filter GetCircles patch.

Algorithms
61

As far as the documentation mentions that computation is done in grid (index) space, it doesn’t seem a problem (to me) to have different spacings.

62

FYI, the helper functions GetCenterPoint() and SetCenterPoint(point) have been added to the master branch, this week:

As well as a notice that GetCenterPoint() yields pixel grid coordinates, when used to retrieve the center of a circle from GetCircles():

A change of behavior can be observed as follows, when the algorithm has detected a circle centered at pixel index [6, 8]:

const itk::SpatialObject<2>* const object = filter->GetCircles().front();
std::cout
  << "To Parent Offset: " << object->GetObjectToParentTransform()->GetOffset()
  << "\nTo World Offset: " << object->GetObjectToWorldTransform()->GetOffset()
  << std::endl;

This produced the following output before the patch:

To Parent Offset: [6, 8]
To World Offset: [0, 0]

Whereas now, after the patch, it produces the following output:

To Parent Offset: [0, 0]
To World Offset: [6, 8]

But of course, instead of calling GetOffset() on whichever transform, to retrieve the center point of a circle, ITKv5 users should just call EllipseSpatialObject::GetCenterPoint() :smiley:

2 Likes
63

@matt.mccormick, @tim-evain I’m glad that you like the patch https://github.com/Kitware/ITK/commit/842f15a86221741f9a52f424907e154d3dbb24bb However, I still don’t fully understand the consequences of having replaced GetObjectToParentTransform()->m_Offset by GetObjectToWorldTransform()->m_Offset as storage for the center coordinates.

Before the patch (for example, with ITK 4.13), it could be checked that the center of a circle is inside its spatial object by calling spatialObject->IsInside(center) as follows:

// Create an image with a circle centered at [6, 8]:
enum { centerX = 6, centerY = 8 };
const auto image = itk::Image<unsigned>::New();
image->SetRegions({ 16, 16 });
image->Allocate(true);
image->SetPixel({ centerX, centerY }, 1);
image->SetPixel({ centerX, centerY - 1 }, 1);
image->SetPixel({ centerX, centerY + 1 }, 1);
image->SetPixel({ centerX - 1, centerY }, 1);
image->SetPixel({ centerX + 1, centerY }, 1);

const auto filter = 
  itk::HoughTransform2DCirclesImageFilter<unsigned, unsigned, double>::New();
filter->SetInput(image);
filter->Update();

// GetCircles() finds a circle of radius 1, centered at [6, 8].
const auto& spatialObject = filter->GetCircles().front();
spatialObject->ComputeObjectToWorldTransform();

const double center[] = { centerX, centerY };
const bool isInside = spatialObject->IsInside(center);

std::cout << (isInside ?
  "OK: The center is inside." :
  "ERROR: The center is not inside!") << std::endl;

Before the patch, this example would print “OK: The center is inside.” Now, after the patch, the example prints “ERROR: The center is not inside!”. Is that OK? Or should IsInside be used in a different way?

Update 26 Feb 2018: The code example is also at https://gist.github.com/N-Dekker/b276d312dc37958b6657f889f2687ba4

1 Like
64

@Niels_Dekker Thanks for looking into this and sharing the example code. :+1:

I looked into it, and after

  this->GetModifiableObjectToWorldTransform()->SetOffset(point - originPoint);

we need to call

  this->ComputeObjectToParentTransform();

so the ObjectToParentTransform is updated after we modified the ObjectToWorldTransform. This is used by IsInside.

Also, we should consider re-naming the method from SetCenterPoint to SetCenter to be consistent with the rest of the toolkit.

65

Thanks, Matt. I did not know about ComputeObjectToParentTransform(). But then, what is the difference between:

Old approach (see also the very first version of GetCircles(), Julien Jomier, 2003):

spatialObject->GetObjectToParentTransform()->SetOffset(center);
spatialObject->ComputeObjectToWorldTransform();

New approach (after the patch):

spatialObject->GetModifiableObjectToWorldTransform()->SetOffset(center);
spatialObject->ComputeObjectToParentTransform();

Why do you think the new approach is better than the old one? Note that the center is currently still expressed in pixel grid coordinates (as it was before).

66

There are two transformation’s inside itk::SpatialObject, related to ObjectToParent and ObjectToWorld. They are related, but two are kept for performance reasons. ComputeObjectToWorldTransform needs to be called in EllipseSpatialObject<TDimension>::GetCenterPoint() so they are properly in sync.

This is a bug, which will cause issues for images with non-unit spacing or non-identity orientation. These lines:

should be replaced by a call to itk::Image::TransformIndexToPhysicalPoint.

It is important to operate in physical space to avoid bugs related to non-unit spacing in images, non-identity orientation of images, images with different sizes, images with non-null origins, comparisons with circles that exist in physical space, meshes that exist in physical space, bounding boxes that exist in physical space, point sets that exist in physical space, etc. ITK’s powerful architecture helps avoid these issues by consistently defining operations and values in physical space. This is why the SpatialObject framework was created.

1 Like
67

Thanks, Matt. If I understand correctly, there are three issues now related to the center of the circles found by HoughTransform2D GetCircles():

  1. ComputeObjectTo… functions should still be called inside either or both GetCenterPoint + SetCenterPoint
  2. HoughTransform2D GetCircles() should include spacial information from the input image.
  3. GetCenterPoint/SetCenterPoint should be renamed to GetCenter/SetCenter.

Do you think each of these issues should be fixed (preferably?) by the first release of ITKv5?

68

:+1: Yes

69

Thanks for merging my proposed patch, which added a ComputeObjectToParentTransform() call to EllipseSpatialObject::SetCenterPoint(point):
https://github.com/Kitware/ITK/commit/b040480ea3f4f86c6f361e57021371c9bcd6f7c1

Do you still think that a call to ComputeObjectToWorldTransform() needs to be added to EllipseSpatialObject::GetCenterPoint() now? If so, can you please show me some test code that currently produces the wrong results, and would be fixed by adding ComputeObjectToWorldTransform() to GetCenterPoint()?

1 Like
70

Thanks for submitting the patch, @Neils_Dekker! That fix was the only place we needed to call ComputeObjectToWorldTransform.

71

You’re welcome, Matt! You mean ComputeObjectToParentTransform(), right? At least, that’s the one I added to SetCenterPoint(point), with ENH: IsInside(point) made easier, especially for HoughTransform circles · Kitware/ITK@210ba16 · GitHub

72

Correct.

1 Like
73

That brings us to the next issue: HoughTransform2D GetCircles() should include spacial information from the input image.

You already suggested calling itk::Image::TransformIndexToPhysicalPoint to convert the found index (2-D) of the center to the physical center point. Would it also be useful (or even necessary) to call itk::SpatialObject::SetSpacing, to copy the pixel spacing from the input image into each found circle? Something like this:

 Circle->SetSpacing(inputImage->GetSpacing());

But then, what would be the actual effect of calling SetSpacing, to the user?

Update, March 5, 2018: The ability to directly pass the spacing of an itk::Image to an itk::SpatialObject was added by Luis Ibanez, 2004-01-24 https://github.com/Kitware/ITK/commit/6956f3390162bbce004704d411ded919a3883fcf but then again removed by Stephen Aylward, 2005-03-09, https://github.com/Kitware/ITK/commit/2b536fbe9d81d74fdf7af42e0950f2a18c5eb861#diff-46707908d0b621a26a7d7943c70bcae1 However, it is still possible to do spatialObject->SetSpacing( inputImage->GetSpacing().GetDataPointer() ).

74

Good questions, @Niels_Dekker

Since the EllipseSpatialObject is defined by its Center and Radius, which are both defined in terms of physical space units, this is all we need to set. We do not need to set the Spacing, whose meaning is ill-defined for an ellipse / circle.

75

OK, thanks for your explanation, @matt.mccormick . Well, it does surprise me somewhat, as in the past, a use case that did actually try to set the spacing of an EllipseSpatialObject was found convincing enough to improve spacing support for SpatialObject. As Luis Ibanez replied to Corinne Mattmann:
https://public.kitware.com/pipermail/insight-users/2004-January/006382.html

Anyway, do you agree that the radii of an object detected by GetCircles() could be as follows:

  const double radius = m_RadiusImage->GetPixel(indexOfMaximum);
  const double radii[] = { radius * spacing[0], radius * spacing[1] };
  Circle->SetRadius(radii);

Instead of this line of code:

Note that it would imply that a circle detected by GetCircles() might have two different radii. Would that be OK to you?

76

No, a circle only has a single radius, which is defined in physical space.

77

Perhaps throw an exception if spacing along axes is different?

78

This is not necessary if the radius is in physical space units.

79

If I remember correctly, processing is in index space, not physical space, hence the suggestion. Otherwise my suggestion is moot.

80

Since the resulting circles are used as an output, the radius should be in physical space.

1 Like