Gradient sometimes zero (which prevents a gradient decent optimizer from decenting)

[romangrothausmann]

I have experienced many cases where I can extract a path with the IterateNeighborhoodOptimizer but where the GradientDescentOptimizer does not start to decent even though the speed function is continuous (based on a binary segmentation), apparently because the initial local gradient is calculated to be zero (even though start- and end-point, speed function etc are kept the same).
I have a vague feeling that this thresholding:

ITKMinimalPathExtraction/include/itkSingleImageCostFunction.hxx

Lines 145 to 152 in 6ae35b9

	// NOTE: The cost function may undefined / unreachable areas
	// (indicated by very large values) which may skew the gradient.
	// To avoid this skewing effect, we reset gradient values larger
	// than a given threshold.
	if ( itk::Math::abs (derivative[i]) > DerivativeThreshold )
	{
	derivative[i] = 0.0;
	}

to a hard coded value of 15.0:

ITKMinimalPathExtraction/include/itkSingleImageCostFunction.hxx

Line 140 in 6ae35b9

constexpr typename DerivativeType::ValueType DerivativeThreshold = 15.0;

could be the culprit.

Another possible reason I could imagine is that the floating point precision might not suffice under some circumstances to calculate a non-zero gradient even if there is no extrema in the cost function at the current optimizer position.

I use this CLI for the testing:
https://github.com/romangrothausmann/ITK-CLIs/blob/bfb1312142d505cacd6770e4d5acc23475290c8f/min-path_seg.cxx

[romangrothausmann]

Where could I set the precision used for the gradient calculation or is it hard-coded somewhere to float or double?

[romangrothausmann]

While testing that #62 works as expected, even setting DerivativeThreshold of 1e9 did not solve the issue reported here, that under some (unknown) circumstances no gradient decent optimizer starts to run, apparently because the computed gradient is zero even though the IterateNeighborhoodOptimizer does find lower values even with the same step-size.
@thewtex Any idea what could cause this? Can you re-open the issue (I can't)?

[romangrothausmann]

It seems float is used by default for the gradient calculation

ITKMinimalPathExtraction/include/itkPhysicalCentralDifferenceImageFunction.h

Lines 41 to 43 in 2ff1ab5

	template <
	class TInputImage,
	class TCoordRep = float >

and to my understanding there is no direct way to override this default when instantiating SpeedFunctionToPathFilter nor by separate instantiating its SingleImageCostFunction, because the PhysicalCentralDifferenceImageFunction

ITKMinimalPathExtraction/include/itkSingleImageCostFunction.h

Line 103 in 2ff1ab5

using GradientImageFunctionType = PhysicalCentralDifferenceImageFunction< ImageType, CoordRepType >;

cannot be re-set.
However, the final default that gets used seems to come from the superclass
CostFunctionTemplate< TInternalComputationValueType >::ParametersValueType

ITKMinimalPathExtraction/include/itkSingleImageCostFunction.h

Line 91 in 2ff1ab5

using CoordRepType = Superclass::ParametersValueType;

Is there any way to set this to double when using SpeedFunctionToPathFilter in a program (e.g. in min-path_seg) without modifying the code of this module by adding another template parameter for this?
Could the limited precision be the source for this issue?

[thewtex]

@romangrothausmann we could also change the default to double if you find that it corrects behavior.

[romangrothausmann]

Many thanks @thewtex for your reply. How can I specify to use double in this case? I.e. how can I specify the template parameter TInternalComputationValueType of the superclass
when it is inherited

ITKMinimalPathExtraction/include/itkSingleImageCostFunction.h

Lines 52 to 53 in 2ff1ab5

	class ITK_EXPORT SingleImageCostFunction :
	public SingleValuedCostFunction

I cannot find where the default is set for this. To my understanding that must happen somewhere, but
there seems no default specification for template< typename TInternalComputationValueType >.
Or what would be other hacky ways that would suffice for testing?

[thewtex]

It looks like the default is already double:

https://github.com/InsightSoftwareConsortium/ITK/blob/99aaeaf1a553866f1e36d4ec363d2a37e1377165/Modules/Numerics/Optimizers/include/itkCostFunction.h#L67

[romangrothausmann]

Many thanks @thewtex for pointing that out. I was not aware of such a way to define a default and was only used to something like

ITKMinimalPathExtraction/include/itkPhysicalCentralDifferenceImageFunction.h

Line 43 in 2ff1ab5

class TCoordRep = float >

@thewtex could you re-open the issue, because it does not seem to be resolved and since the default already is double I've currently no idea what else could be the source for the problem, do you?

[thewtex]

I am not sure -- unfortunately, I do not have time to dig into it at the moment.

[romangrothausmann]

Edited the issue description, now mentioning that the speed image is based on a binary segmentation image (dilated by a parabolic SE in a simple case), so noise should not be the issue.

[thewtex]

How about using an inverse distance map generated from the segmentation for the speed image?

[romangrothausmann]

Thanks @thewtex for that idea, which would be much quicker to generate than my current approach using a fast-marching map of the 3D skeleton (inside the segmentation) squeezed into the range [0; 1] by a sigmoid. However, this has the advantage that the speed along the center line is always close to 1 independent of the local extent of vessels (which vary by a factor of 1000 in our dataset).

Still, I think that another speed map would only represent a work-around for the issue here, because I would expect the local gradient not to equal zero if IterateNeighborhoodOptimizer does start to propagate at the start-point, even reaching the end-point.

[romangrothausmann]

Just to mention that I use this CLI for the testing:
https://github.com/romangrothausmann/ITK-CLIs/blob/bfb1312142d505cacd6770e4d5acc23475290c8f/min-path_seg.cxx