Hough Transform 2D Circles Image Filter GetCircles patch.

tim-evain · February 5, 2018, 9:46am

@Niels_Dekker
I’m not on Gerrit, but your patch seems fine to me.

Niels_Dekker · February 16, 2018, 1:51pm

@tim-evain So far, the patch ENH: Added HoughTransform2D Circle class, simplified GetCircles() return type has not been accepted at http://review.source.kitware.com/#/c/23139

Now I just submitted a new patch, ENH: Added GetCenterPoint and SetCenterPoint to EllipseSpatialObject which would also simplify using GetCircles(), a little bit, at least! Please have a look: http://review.source.kitware.com/#/c/23188

I wouldn’t say that the new patch entirely supersedes the old one, but I hope it can be accepted and merged more easily.

tim-evain · February 17, 2018, 3:32pm

Hi Niels,

My 2 cents:

As Matt said, “auto” shouldn’t be used, there are still some occurences to change to explicit type (lines 55 / 62 / 79 in itkEllipseSpatialObject.hxx for example).
Following the idea that the spatial object is referred in world coordinates, the Hough transform GetCircles() is not consistent as calculated center/radius are in (continuous) indexes.
Edit : there is also a build error on itkEllipseSpatialObject.hxx line 89 “transform->SetOffset(point - originPoint)” due to constness I think.

Tim

Niels_Dekker · February 17, 2018, 5:03pm

Thanks for your feedback, Tim-Evain.

Does ITK ban the C++11 auto keyword? I hope not.

What do you mean by “continuous indexes”? The coordinates of the centers of the circles produced by GetCircles() are always whole numbers, even while they are stored as floating point numbers, in ellipseSpatialObject->GetObjectToParentTransform()->m_Offset. This has always been the behavior of GetCircles(), from the very first commit (about 15 years ago).

Good catch! I’ll try to fix it with my next patch set.

The biggest problem seems to me that replacing the call to GetObjectToParentTransform() by GetObjectToWorldTransform(), as I did with Patch Set 3, would silently break user code. What do you think?

tim-evain · February 17, 2018, 8:28pm

To be honest I don’t know what are the official guidelines on this. But as much as “auto” is a great feature for prototyping and testing, it goes fundamentally against the C(++) spirit of strong typing. In open-source project like ITK, it hinders maintenance for contributors and understanding for any users. Going with explicit type also ensure that the code is really doing what the programmer expects it to do: you are worried about silently breaking compatibility. That’s exactly what could happen without careful use of “auto”.

On the change of ObjectToParent to ObjectToWorld, I let the officials take the decision.

Center is a standard index. Radius are continuous indexes. My point is that they are stored as is into the spatial object, but users are gonna expect world ones:

/* Returns the center point (in world coordinates) of the ellipse */
  PointType GetCenterPoint()  const;

Niels_Dekker · February 18, 2018, 2:47pm

Thanks for your reply, @tim-evain
The discussion on auto continues at Using the C++11 auto keyword in ITKv5 code

Almost everything in HoughTransform2DCircles is calculated and expressed in pixel coordinates, not in world coordinates. If we want to keep it that way, a SpatialObject may not be the right data structure to store the result of GetCircles().

However, if we would like to have all input parameters and output of HoughTransform2DCircles to be expressed in world coordinates, should the filter still accept an input image for which spacing[0] != spacing[1]?

tim-evain · February 19, 2018, 10:08am

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 · February 23, 2018, 1:15pm

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

Niels_Dekker · February 24, 2018, 1:16pm

@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

matt.mccormick · February 26, 2018, 7:17pm

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

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.

Niels_Dekker · February 27, 2018, 12:43pm

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

matt.mccormick · February 28, 2018, 5:16pm

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:

github.com

InsightSoftwareConsortium/ITK/blob/e41790229b92fd98865037153ac631e68ac0e9fa/Modules/Filtering/ImageFeature/include/itkHoughTransform2DCirclesImageFilter.hxx#L251-L253


      
          CircleType::PointType center;
          center[0] = indexOfMaximum[0];
          center[1] = indexOfMaximum[1];

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 · February 28, 2018, 8:47pm

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

ComputeObjectTo… functions should still be called inside either or both GetCenterPoint + SetCenterPoint
HoughTransform2D GetCircles() should include spacial information from the input image.
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 · February 28, 2018, 9:26pm

Yes

Niels_Dekker · March 2, 2018, 2:28pm

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

matt.mccormick · March 2, 2018, 2:45pm

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

Niels_Dekker · March 2, 2018, 2:51pm

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

matt.mccormick · March 2, 2018, 3:48pm

Correct.

Niels_Dekker · March 2, 2018, 4:34pm

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

matt.mccormick · March 5, 2018, 4:55pm

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.