Hough Transform 2D Circles Image Filter GetCircles patch.

dzenanz · December 23, 2017, 6:48pm

It is good that you noticed the error. It would be ideal if you also fixed it by submitting a patch

Niels_Dekker · December 29, 2017, 3:54pm

OK, here it is: BUG: GaussianDerivativeImageFunction should use image spacing consistently http://review.source.kitware.com/#/c/22954/

The fix is small (just moving a single line of code), but I think it’s a serious bug, that’s why I raised attention for this issue, here at the discussion forum.

Niels_Dekker · December 29, 2017, 7:51pm

For the record, the following little example shows how increasing the spacing of the input image might affect the result of GaussianDerivativeImageFunction, now that my bug fix on this subject https://itk.org/gitweb?p=ITK.git;a=commit;h=4119f2f1 has been merged:

#include "itkGaussianDerivativeImageFunction.h"
#include <iostream>

int main()
{
  for (double spacing = 0.25; spacing < 5; spacing *= 2)
  {
    typedef itk::Image<unsigned int> ImageType;
    const auto image = ImageType::New();
    image->SetRegions({ 2, 1 });
    image->Allocate(true);
    image->GetBufferPointer()[0] = 1;
    image->SetSpacing(spacing);

    typedef itk::GaussianDerivativeImageFunction<ImageType> ImageFunctionType;
    const auto imageFunction = itk::GaussianDerivativeImageFunction<ImageType>::New();
    imageFunction->SetInputImage(image);
    const auto result = imageFunction->EvaluateAtIndex({ 0, 0 });

    std::cout << "Spacing = " << spacing << "; result = " << result << std::endl;
  }
}

Output (from Visual C++ 2017 64-bit, on Win10):

Spacing = 0.25; result = [0.484617, 0]
Spacing = 0.5; result = [0.882497, 0]
Spacing = 1; result = [1.21306, 0]
Spacing = 2; result = [0.541341, 0]
Spacing = 4; result = [0.0026837, 0]

So in this case, the result goes up and down

But I’m not sure if these values should be considered the “ground truth”, as there appears to be another major GaussianDerivativeImageFunction bug, which I don’t know how to fix properly (WIP: Removed blurring from GaussianDerivativeImageFunction::EvaluateAtIndex, http://review.source.kitware.com/#/c/22926).

blowekamp · December 29, 2017, 8:15pm

Your test case is not a good one. It is a 1x2 image where you are taking a gaussian derivative. It is literally a corner case. This involves the boundary condition which may show other problems.

I’d suggest something like a 65x65 image and to look at the center pixel.

Regarding not “trusting” this implementation of the function, it is only used in ITK proper by the Hough. When I implemented a 3D Hough, I still used this functions and I got good results.

If you are looking for an alternative, in review there is: https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Nonunit/Review/include/itkDiscreteGaussianDerivativeImageFunction.h

I recall looking very close over that function sometime ago, so I trust it.

Niels_Dekker · December 29, 2017, 9:38pm

@blowekamp OK, I just tried a 65x65 image. However, the following does produce the very same output as with my original 2x1 image:

for (double spacing = 0.25; spacing < 5; spacing *= 2)
{
  typedef itk::Image<unsigned int> ImageType;
  const auto image = ImageType::New();
  enum { imageSize = 65 };
  image->SetRegions({ imageSize, imageSize });
  image->Allocate(true);
  image->SetPixel({ imageSize / 2, imageSize / 2 }, 1);
  image->SetSpacing(spacing);

  typedef itk::GaussianDerivativeImageFunction<ImageType> ImageFunctionType;
  const auto imageFunction = itk::GaussianDerivativeImageFunction<ImageType>::New();
  imageFunction->SetInputImage(image);

  const auto result = imageFunction->EvaluateAtIndex({ imageSize / 2 + 1, imageSize / 2 });
  std::cout << "Spacing = " << spacing << "; result = " << result << std::endl;
}

Note that with default sigma (1.0) and default extents (1.0), the kernel radius is only 1, so I think a smaller test image would be fine as well anyway… For the record, the kernel size is set here, in RecomputeGaussianKernel(): https://github.com/Kitware/ITK/blob/v4.13.0/Modules/Core/ImageFunction/include/itkGaussianDerivativeImageFunction.hxx#L169

If you think that the current result of itk::GaussianDerivativeImageFunction is still OK for Hough, we could make it twice as fast (at least) by removing the entire blurring step (as in http://review.source.kitware.com/#/c/22926). What do you think?

The alternative would be to replace GaussianDerivativeImageFunction by something else, indeed. Thanks for suggesting itkDiscreteGaussianDerivativeImageFunction. I’ll have a look at it later.

Niels_Dekker · January 2, 2018, 9:32pm

FYI, it appears that the patch I just submitted PERF: NeighborhoodOperatorImageFunction avoids copy ConstNeighborhoodIterator http://review.source.kitware.com/#/c/22959/ significantly improves the performance of the HoughTransform2DCircles filter. I reran the test that I posted at Dec 10, 2017 at Hough Transform 2D Circles Image Filter GetCircles patch. and saw a performance improvement of more than 30% (duration going from 45 seconds before the patch, to less than 30 seconds).

tim-evain · January 24, 2018, 10:38am

Hi all,

I bring back this thread as it is closely related to my “bug” discovery, but if an admin feels like it should be a new post, feel free to do it.

I’ve noticed that GetCircles method is quite tricky about the center of found circles.
When calling:

  CircleType::VectorType center;
  center[0] = indexOfMaximum[0];
  center[1] = indexOfMaximum[1];
  Circle->GetObjectToParentTransform()->SetOffset(center);
  Circle->ComputeBoundingBox();

the GetObjectToParentTransform method of SpatialObject:

template< unsigned int TDimension >
typename SpatialObject< TDimension >::TransformType *
SpatialObject< TDimension >
::GetObjectToParentTransform(void)
{
  return static_cast< TreeNodeType * >(
           m_TreeNode.GetPointer() )->GetNodeToParentNodeTransform();
  //return m_ObjectToNodeTransform.GetPointer();
}

return a pointer to the NodeToParentNodeTransform in place of the ObjectToParentTransform, which is quite confusing when trying to get the center back. Is that intended ?

Tim

Niels_Dekker · January 24, 2018, 11:05am

Thanks Tim. I do think HoughTransform2DCirclesImageFilter::CircleType is overly complicated. It appears that the proper way to retrieve the center from a CircleType object is something like circle.GetObjectToParentTransform()->GetOffset(). I would like it to be easier.

Anything as simple as the following would be just fine to me:

struct CircleType
{
  double radius;
  double center[2];
};

Of course, such a change would break old ITK user code, but maybe it is acceptable for ITK 5? What do you think?

tim-evain · January 24, 2018, 11:45am

Even if it works, it is error prone as you are getting in fact the NodeToParentNodeTransform, not the ObjectToParentTransform…

That would be simpler, but I can see the interest of using an EllipseSpatialObject, in case of concentric circles. I would advocate for a fixed processing with SpatialObject, then we can access circle information by its bounding box center/bounds.

Tim

Niels_Dekker · January 24, 2018, 1:17pm

Do you need to use a circle, created by ITK’s Hough Transform, as EllipseSpatialObject, at the moment? I don’t, I just need to have its center and its radius.

For those cases where the circle would really need to be represented as an EllipseSpatialObject, a function could be provided (outside of ITK’s Hough Transform), something like CreateEllipseSpatialObjectForCircle(double radius, const double center[2]), returning an EllipseSpatialObject. Right?

Of course it might still be useful to improve itk::EllipseSpatialObject, if that’s what you’re suggesting. I think that could be done independently of the question how to improve the return type of itk::HoughTransform2DCirclesImageFilter::GetCircles().

tim-evain · January 24, 2018, 2:20pm

For my personal usage no.

If you provide a custom structure into the Hough filter, I think it should ship the conversion function, but it could be out of the core processing for sure.

Niels_Dekker · February 1, 2018, 10:02pm

@tim-evain Can you please have a look at my work in progress, WIP: Added Circle class to simplify HoughTransform2D GetCircles() at http://review.source.kitware.com/#/c/23139/ ?

What do you think? As you see, it’s only simplifying the circle representation, not the itk::EllipseSpatialObject template. If you still think EllipseSpatialObject should be improved as well, I would suggest you to write another proposed patch!

tim-evain · February 2, 2018, 3:20pm

Hi Niels!

I think this is a great improvement over the previous system .
Some thoughts:

On “ToDo 3” : maybe start with the definition of Circle into the Hough file, and see afterward if there is really a need for other filters to use it ?
Circles are found by an index center, but defined as a point afterward, leading to a confusing representation. An user that doesn’t check the code will think that the circle is defined in world coordinates, but the radius and center are in fact in continuous indexes. I think it should be harmonized, either by making the center an index, or putting all values in world coordinates (I think the latter is preferable).
On “ToDo 2” : it depends on the previous choice. If the circle is in index system, there is no need for template.
Since the return type of GetCircles() will change, I wonder if a vector is not preferable over a list. The syntax to access circles with lists is heavier in my opinion.

Anyway thanks for the time you put into making this!

Niels_Dekker · February 2, 2018, 8:17pm

@tim-evain Please check the new Patch Set that I just submitted:
ENH: Added HoughTransform2D Circle class, simplified GetCircles() return type
http://review.source.kitware.com/#/c/23139/4

It should address most of your remarks:

It placed the Circle class nested inside the HoughTransform2DCircles class
It represented the center of the circle as an Index<2>, to clarify that the (x,y) are not yet world coordinates. (Switching to floating point world coordinates might be a future improvement. Although for my personal use cases, the Index<2> is just appropriate.)
It switched from std::list to std::vector. I wanted to do so before, but I needed somebody to encourage me

Note that the current HoughTransform2DCircles interface allows the user to specify the type of the pixels of the radius image, as a template argument, TRadiusPixelType. It seemed reasonable to me to use this as the type of the radius of the Circle class as well.

If you have access to http://review.source.kitware.com please review my patch directly at http://review.source.kitware.com/#/c/23139

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]?