smoothstep returning out of range values for 0 length steps

[Skwint]

Godot version

4.0.dev (merge #68060)

System information

Windows 11

Issue description

In math_funcs.h the two smoothstep functions return a value clamped between 0.0 and 1.0 except when the from and to values are almost equal, in which case they return the value of from, which could be 7 million for all we know. I'm confident they should return a value between 0.0 and 1.0 instead. Personally I think 0.5 would be a decent compromise.

Steps to reproduce

I'm not reproducing this I'm just eyeballing the code. It's only 3 lines long.

Minimal reproduction project

No response

[Lielay9]

It would be nice if this was also reflected on the docs. I don't think it's currently 100% up to date with how the function works since the to can also be lower than from.

Returns the result of smoothly interpolating the value of x between 0 and 1, based on the where x lies with respect to the edges from and to.
The return value is 0 if x <= from, and 1 if x >= to. If x lies between from and to, the returned value follows an S-shaped curve that maps x between 0 and 1.
...

(There's a typo too btw: ...based on the where x... should be ...based on where the x...)

From a glance, it seems that at least in OpenGL from is expected to be lower.

Results are undefined if edge0 ≥ edge1.

So if we expect the same (which the current doc alludes to), the function could fall back to conditionals when the edges are close to equal. Something like

if (is_equal_approx(p_from, p_to)) {
	if (p_s <= p_from)
		return 0; 
	if (p_s >= p_to)
		return 1;
	return 0.5;
}

Not sure if it's any better though.

[Skwint]

I think I agree with the function you wrote, except the 0.5 result should almost never happen so you can drop an if statement and always give 0 or 1.

[clayjohn]

I guess the check for equality is to avoid a division by zero in the next line:

godot/core/math/math_funcs.h

Lines 392 to 405 in 604abb4

	static _ALWAYS_INLINE_ double smoothstep(double p_from, double p_to, double p_s) {
	if (is_equal_approx(p_from, p_to)) {
	return p_from;
	}
	double s = CLAMP((p_s - p_from) / (p_to - p_from), 0.0, 1.0);
	return s * s * (3.0 - 2.0 * s);
	}
	static _ALWAYS_INLINE_ float smoothstep(float p_from, float p_to, float p_s) {
	if (is_equal_approx(p_from, p_to)) {
	return p_from;
	}
	float s = CLAMP((p_s - p_from) / (p_to - p_from), 0.0f, 1.0f);
	return s * s * (3.0f - 2.0f * s);
	}

I would probably just clamp p_to - p_from to a small positive value like:

float range = MAX(p_to - p_from, CMP_EPSILON);
float s = CLAMP((p_s - p_from) / range, 0.0f, 1.0f); 
return s * s * (3.0f - 2.0f * s); 

[Skwint]

That's going to give you a headache if p_to is less than p_from ...

[Lielay9]

That's going to give you a headache if p_to is less than p_from ...

That could be changed with one extra if. Taking inspiration from the latter function, move_toward:

...
float range = abs(p_to - p_from) > CMP_EPSILON? p_to - p_from : SIGN(p_to - p_from) * CMP_EPSILON;
...

Admittedly, it's getting quite messy...

[Lielay9]

To stir the pot a little more, it looks like the shaders smoothstep defaults to 0 (at least on my machine).

COLOR.rgb = vec3(smoothstep(1.0, 0.9999994, 0.0)) // white color rect
COLOR.rgb = vec3(smoothstep(1.0, 0.9999995, 0.0)) // black color rect

[Skwint]

If the values are 'equal' we have no reliable way of knowing whether they are ascending or descending (if they do differ it is likely because the 'same' value was arrived at via a different sequence of equations and the delta is nothing more than rounding errors) so we don't know whether greater or less than should be the 0.
I think it really depends on the use case.
If the result is supposed to fade a colour (as a simple example) then a delta of 0 means an instant colour change with no fade, but what we will get is no colour anywhere (0), all colour everywhere (1), half colour everywhere (0.5) or an instant colour change that is in the wrong direction 50% of the time and likely to cause flicker (if statements).

A lot of times what I really want here is an error, at least in debug mode, because it means I've screwed up. I should not be calling this function with these values.

Whatever is chosen, however, if I multiply it by 255 and insert it into a colour channel it should not be possible for it to overflow into a different channel. If I mutiply the return value by the length of an array it should always give a valid index. I've made a mistake, but I'd like to fail gracefully.

[clayjohn]

That's going to give you a headache if p_to is less than p_from ...

That's not a valid input. In such a case we should just error and return 0 I don't think it is a good idea to try to gracefully handle invalid inputs. But thanks for pointing that out. We should indeed be handling this case and not silently returning garbage

To stir the pot a little more, it looks like the shaders smoothstep defaults to 0 (at least on my machine).

COLOR.rgb = vec3(smoothstep(1.0, 0.9999994, 0.0)) // white color rect
COLOR.rgb = vec3(smoothstep(1.0, 0.9999995, 0.0)) // black color rect

For shaders it is undefined behaviour, see: https://registry.khronos.org/OpenGL-Refpages/gl4/html/smoothstep.xhtml In my experience the result can really be anything. Some devices handle it gracefully while others give you numbers that make no sense.

For us, I would just add an error condition when p_from > p_to

Final code for float version would be:

static _ALWAYS_INLINE_ float smoothstep(float p_from, float p_to, float p_s) { 
    ERR_COND_V(pfrom >= p_to, 0.0);

    float range = MAX(p_to - p_from, CMP_EPSILON);
    float s = CLAMP((p_s - p_from) / range, 0.0f, 1.0f); 
    return s * s * (3.0f - 2.0f * s); 
 } 

[Lielay9]

For shaders it is undefined behaviour, see: https://registry.khronos.org/OpenGL-Refpages/gl4/html/smoothstep.xhtml In my experience the result can really be anything. Some devices handle it gracefully while others give you numbers that make no sense.

For us, I would just add an error condition when p_from > p_to

Final code for float version would be:

static _ALWAYS_INLINE_ float smoothstep(float p_from, float p_to, float p_s) { 
    ERR_COND_V(pfrom >= p_to, 0.0);

    float range = MAX(p_to - p_from, CMP_EPSILON);
    float s = CLAMP((p_s - p_from) / range, 0.0f, 1.0f); 
    return s * s * (3.0f - 2.0f * s); 
 } 

If my memory serves correctly there was some discussion somewhere about the doc being a bit unspecific about the undefined behaviour of smoothstep. The only real error case is when from and to are equal or close to equal. When to is less than from the function just "flips" for lack of a better word. I've seen this "trick" used quite often in shaders. E.g:

smoothstep(max, min, t)
// instead of
1.0 - smoothstep(min, max, t)

That said if Wikipedia is to trust, the definition does assume from being the smaller of the two and at least, for example, Unity disallows using lesser to (defaults to returning 0). Personally, I favour the "flipping" behaviour and would only error out on equal inputs, but it's not really that big of a deal either way. It might even be better to follow the norm if for some reason the shader implementation changes in the future.

[thygrrr]

I need to chime in here (how has this been untouched for a year?)

The function is also quite incorrectly documented.

SmoothStep doesn't "interpolate a value", it returns a value from 0 to 1 (except in the edge case discussed here). I like the suggestion of using 0.5 if the values are approximately equal. If the extra cases of the value being inside and outside the range can be added (so logical instead of analytical solution), that would be good instead of trying to make the code maximally branchless like for a GL.

The documentaion error is also there in the C# version:

    /// <summary>
    /// Returns a number smoothly interpolated between <paramref name="from" /> and <paramref name="to" />,
    /// based on the <paramref name="weight" />. Similar to <see cref="M:Godot.Mathf.Lerp(System.Single,System.Single,System.Single)" />,
    /// but interpolates faster at the beginning and slower at the end.
    /// </summary>
    /// <param name="from">The start value for interpolation.</param>
    /// <param name="to">The destination value for interpolation.</param>
    /// <param name="weight">A value representing the amount of interpolation.</param>
    /// <returns>The resulting value of the interpolation.</returns>
    public static float SmoothStep(float from, float to, float weight)
    {
      if (Mathf.IsEqualApprox(from, to))
        return from;
      float num = Math.Clamp((float) (((double) weight - (double) from) / ((double) to - (double) from)), 0.0f, 1f);
      return (float) ((double) num * (double) num * (3.0 - 2.0 * (double) num));
    }

[Malcolmnixon]

This function is definitely wrong in it's current implementation; but a fix has been held up for quite a while trying to find a good implementation. My thoughts are as follows:

Degenerate Value (when from == to)

It almost certainly makes no difference whether we return 0.0, 0.5, or 1.0, as floats aren't perfectly accurate except in extremely rare/contrived cases. Practically speaking no-one will know how far x is in the zero-width range anyway. The only thing we can't do is return a value outside of the 0.0 - 1.0 range this function is supposed to generate.

If I had to pick an answer then I'd pick 0.5 - as it's both splitting the difference, and the solution to the somewhat similar Grandi's series.

Range Order

The current implementation supports both from < to and from > to. We can argue about whether it should have, or whether OpenGL does it a different way; however Godot has been shipping with smoothstep working with positive and negative ranges (just like inverse_lerp and remap), and it's too late to change it now without risking breaking existing users.

Proposed Implementation

I propose the following implementation based on the thoughts above:

	static _ALWAYS_INLINE_ double smoothstep(double p_from, double p_to, double p_s) {
		if (is_equal_approx(p_from, p_to)) {
			return 0.5;
		}
		double s = CLAMP((p_s - p_from) / (p_to - p_from), 0.0, 1.0);
		return s * s * (3.0 - 2.0 * s);
	}
	static _ALWAYS_INLINE_ float smoothstep(float p_from, float p_to, float p_s) {
		if (is_equal_approx(p_from, p_to)) {
			return 0.5f;
		}
		float s = CLAMP((p_s - p_from) / (p_to - p_from), 0.0f, 1.0f);
		return s * s * (3.0f - 2.0f * s);
	}

Proposed Documentation

Returns a smooth cubic Hermite interpolation between 0 and 1.

When from < to the return value is 0 if x <= from, and '1' if x >= to. if x lies between from and to, the return value follows an S-shaped curve that smoothly transitions from '0' to '1'.

When from > to the function is mirrored and returns '1' if 'x <= to' and '0' if 'x >= from'.

This S-shaped curve is the cubic Hermite interpolator, given by f(y) = 3*y^2 - 2*y^3 where y = (x-from) / (to-from).

[BastiaanOlij]

if (is_equal_approx(p_from, p_to)) {
	return 0.5;
}

It actually wrong, this means the output is always 0.5 regardless of p_s, but the failure situation is only if p_s lies between from and to. If p_s lies outside of the range, it should be either 0.0 or 1.0.

So correct would be (and i'm ignoring for a minute Malcolms suggestion to reverse the curve if p_to < p_from):

if (is_equal_approx(p_from, p_to)) {
	return p_s < p_from ? 0.0 : 1.0;
}

I really don't think there is a point of returning 0.5 if p_s == p_from, but if that is important it can be expanded to:

if (is_equal_approx(p_from, p_to)) {
	if (p_s < p_from) {
		return 0.0;
	} else if (p_s > p_to) {
		return 1.0;
	} else {
		return 0.5;
	}
}

[Malcolmnixon]

Ah yes, in the degenerate case you can consider EVERY value to be outside the range, so it's got to be 0.0 or 1.0.

If we treat the zero-sized range as "normal" polarity (not inverted) then the correct code is as you stated:

if (is_equal_approx(p_from, p_to)) {
	return p_s < p_from ? 0.0 : 1.0;
}