Fast parallel lines and vectors test

Picture showin the similiar ratio for 2 triangles.

Yesterday I decided to stop being lazy in math, as I was recently invited to become part of the Irrlicht team I started with small tasks regarding math utils.

However, I basically copy-pasted few algorithms (of course by linking original author) and readapted some methods to be more specific so that users can choose functions more suitable to their uses cases. It was a small contribution, but a lazy one!

I stopped being lazy when I had to implement a function for detecting if 2 lines/vectors were parallel. I did not found a ready implementation in C++ that was easy to understand to everyone (some implementations even involved ranking a matrix, which requires some calculus background).

So I finally decided to stop being lazy and to do something as easy as possible to understand and fast at the same time.

Just note that parallel lines have X/Y components that are in the same ratio, like similar triangles.

Picture showin the similiar ratio for 2 triangles.
The catheti in the two triangles A and B have the same ratio and so the same inclination for the hypotenuse
  • If 2 lines are parallel, then the catheti of any triangle subtended by the line on the X axis have the same ratio.
  • Of course, there are some special cases we need to handle, in particular for vertical and horizontal axes, but we can just check if the X or Y components are both 0.

The first version of the method was: not very clever!


template <class T>
class vector2d
{
    bool nearlyParallel ( const vector2d<T> & other ) const
    {
        return
            ( X == 0 && other.X == 0)
        ||  ( Y == 0 && other.Y == 0)
        ||  ( X/Y == other.X / other.Y);
    }
};

So far so good, but in reality, the above code has some serious problems:

  1. It may incur in division by zero
  2. Division is slow
  3. It has to be valid also for integer overload of the class (all Irrlicht’s math have to be valid for most native integer and floating point types)
  4. It is hard to find 2 perfectly parallel vectors with floating point math, in particular, the above method return “false” even in the case where one vector is compared with its normalized version.

You can solve all these problems by using a simple trick:

  • If  X/Y = Z/W, then X*W = Z*Y
  • use a comparison function with some kind of tolerance

template <class T>
class vector2d
{
    bool nearlyParallel ( const vector2d<T> & other ) const
    {
        return
            ( X == 0 && other.X == 0)
            ||  ( Y == 0 && other.Y == 0)
            || equals( X*other.Y, other.X *Y); // comparison and multiplication here
    }
};

The above code could almost work, we have still two issues:

  • We have a boolean expression that is or’ing three different equations, to test parallelity we are forced to execute all the three expressions.
  • The are some cases where we get fake “positives”.

If two vectors are orthogonal and aligned to X, Y axes we get a fake positive. Also, zero vectors cause some troubles, for simplicity we want zero vectors to never be parallel to another vector.

In this verions I fix the issues, allowing for a short circuit in case first expression is false, and by removing fake positives:


template <class T>
class vector2d
{
    bool nearlyParallel( const vector2d<T> & other) const
    {
        return
            equals( X*other.Y, other.X* Y) &&
            (X*other.X + Y*other.Y) != 0;
    }
};

We just have to make sure that only X or only Y are 0 at the same time, if not then the second part of the expression return false.

However, we still have one big problem, that went unnoticed because it is in the implementation of the “equals” function in Irrlicht: it uses an absolute error.


template <>
inline bool equals(const f32 a, const f32 b, const f32 tolerance = ROUNDING_ERROR_f32)
{
    return (a + tolerance >= b) && (a - tolerance <= b);
}

I had then to implement comparison function using the relative error. But here again, the problem with a relative comparison function is that it has some branching and a division:


template <class T>
inline bool equalsRelative( const T a, const T b, const T factor)
{
    if( a==b)
        return true;

    T c = (T) fabs(a);
    T d = (T) fabs(b);

    T relativeError = 0;
    if(c>d)
        relativeError = (T)fabs((c-d)/c);
    else
        relativeError = (T)fabs((c-d)/d);

    return relativeError <= factor;
}

Can we do better? Of course. Ternay operators helps us to remove some branching:


template <class T>
inline bool equalsRelative( const T a, const T b, const T factor = relativeErrorFactor<T>())
{
    T maximum = max_( a, b);
    T minimum = min_( a, b);

    if(maximum == 0)
        return maximum == minimum;
    return  (T) fabs( (maximum-minimum)/maximum) < factor;
}

We still have a division and we cannot multiply by “maximum”  (to get rid of the division) because it can change the sign when taken outside of the absolute value function thus it would require some other changes. Also, we introduced a branch for preventing division by zero.

If we think about why we need relative error, it is needed to determine if the error is going to change effectively our number, basically, if a number is small enough compared to another number, it will not change it when summed.

So I changed approach and I checked if the difference between the numbers was going to be noticed.


template <class T>
inline bool equalsRelative( const T a, const T b, const T factor = relativeErrorFactor<T>())
{
    T maximum = max_( a, b);
    T minimum = min_( a, b); //minimum could be greater in magnitude if it is negative...
    T delta = maximum - minimum;

    return  (factor*maximum + delta == factor*maximum) &&
            (factor*minimum - delta == factor*minimum); // ..so we need to check also "minimum"
}

It looks we are near to an optimal version of the function: if we expand the “delta” in fact


    return  (factor*maximum + (maximum - minimum) == factor*maximum) &&
            (factor*minimum - (maximum - minimum) == factor*minimum);

We can then rearrange variables in order to exploit 4 MAD (MultiplyAndAdd) instructions of the CPU,


    return  (factor*maximum + maximum == factor*maximum + minimum) &&
            (factor*minimum + minimum == factor*minimum + maximum);

At this point, we observe that:

  • Both expressions will be equals only if delta is smaller than a factor of one of the two inputs.
  • One input is always greater than or equal to the other input in magnitude.
  • We can just keep the expression with greater input (if we first find the greater input), but it is hard to tell what is faster without a profiler, and anyone we are already pretty fast.

So here’s the recap of the improvements ported to Irrlicht:

//! returns if a equals b, taking relative error in form of factor
//! this particular function does not involve any division.
template <class T>
inline bool equalsRelative( const T a, const T b, const T factor = relativeErrorFactor<T>())
{
T maximum = max_( a, b);
T minimum = min_( a, b);

return    (maximum*factor + maximum) == (maximum*factor + minimum)
&&(minimum*factor + maximum) == (minimum*factor + minimum); // MAD Wise
}

template <class T>
class vector2d
{
//...

//! check if this vector is parallel to another vector
bool nearlyParallel( const vector2d<T> & other) const
{
return  equalsRelative( X*other.Y, other.X* Y)
&& (X*other.X + Y*other.Y) != 0;
}

//...
};

 

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s