Hough Transform 2D Circles Image Filter GetCircles patch.

(Tim Evain) #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.

(Niels Dekker) #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();
  << "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:

(Niels Dekker) #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->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();

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

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

(Matt McCormick) #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


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.

(Niels Dekker) #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):


New approach (after the patch):


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).

(Matt McCormick) #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.

(Niels Dekker) #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?

(Matt McCormick) #68

:+1: Yes

(Niels Dekker) #69

Thanks for merging my proposed patch, which added a ComputeObjectToParentTransform() call to EllipseSpatialObject::SetCenterPoint(point):

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()?

(Matt McCormick) #70

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

(Niels Dekker) #71

You’re welcome, Matt! You mean ComputeObjectToParentTransform(), right? At least, that’s the one I added to SetCenterPoint(point), with https://github.com/Kitware/ITK/commit/210ba16c02641dfc14d7b1381aa19545fbdb0338

(Matt McCormick) #72


(Niels Dekker) #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:


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() ).

(Matt McCormick) #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.

(Niels Dekker) #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:

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] };

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?

(Matt McCormick) #76

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

(Dženan Zukić) #77

Perhaps throw an exception if spacing along axes is different?

(Matt McCormick) #78

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

(Dženan Zukić) #79

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

(Matt McCormick) #80

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