Explain this line segment intersection code

I've been reading this article on pathfinding over polygon meshes using line-of-sight, and in the article are several code snippets. One of them, the one I'm trying to understand, is a simple function that determines whether two line segments intersect.

Here is that function from the article (in C#, though that isn't important to me):

public static bool LineSegmentsCross(Vector2 a, Vector2 b, Vector2 c, Vector2 d)
{
float denominator = ((b.X - a.X) * (d.Y - c.Y)) - ((b.Y - a.Y) * (d.X - c.X));

if (denominator == 0)
{
return false;
}

float numerator1 = ((a.Y - c.Y) * (d.X - c.X)) - ((a.X - c.X) * (d.Y - c.Y));

float numerator2 = ((a.Y - c.Y) * (b.X - a.X)) - ((a.X - c.X) * (b.Y - a.Y));

if (numerator1 == 0 || numerator2 == 0)
{
return false;
}

float r = numerator1 / denominator;
float s = numerator2 / denominator;

return (r > 0 && r < 1) && (s > 0 && s < 1);
}


Without any explanation, I assume that the function accepts two line segments, a-b and c-d. I know the first part checks whether the two lines are parallel, in which case they don't intersect. After that, I'm not sure. I suspect the next bit checks whether the two segments are out of range of each other? The code talks about numerators and denominators, so maybe it's dealing with slope comparisons?

In any case, I'd like to understand the math used here for determining line intersections. I know there's an algorithm that this code is following, but I don't know what it is. How does the formula used in the first part determine that the lines are parallel? I know it works, but why does it work? And what does the function do after that? Why?

I'd much rather understand how the code works, and then write my own code from that knowledge, than just blindly copy someone else's code.

First, the value called denominator is something that you might call the "two-dimensional cross product" of the vectors b - a and d - c. You can see that it corresponds to what would be the z-component of the cross product if these were 3D vectors. The "2D cross product", like the 3D one, represents the (signed) area of the parallelogram spanned by the two vectors. If the two vectors are parallel or antiparallel, the parallelogram squashes down to nothing and its area is zero.
The two numerators are another pair of 2D cross products. In particular, numerator1 is the cross of the vectors d - c and a - c, that is a vector along one line with a vector pointing from that line toward the other line. So this calculates a parallelogram that spans between the two lines.
At this point it's helpful to sit down and draw a few diagrams, but the key thing to realize is that when the two line segments are barely touching—i.e. when they touch only at an endpoint—the parallelogram spanned by d - c and a - c is either squashed down to nothing, or it's equal to the previous parallelogram (the one calculated in denominator). In other words, in these barely-touching cases, numerator1 is either zero or equal to denominator.
It's not too hard to convince yourself, then, that intermediate cases where the line segments intersect somewhere in the middle result in parallelogram areas somewhere between 0 and denominator, or in other words, numerator1 / denominator should be between 0 and 1.
Then numerator2 is calculated the same way, but just swapping segments. It's necessary to check both directions because there are cases where the line segments don't intersect, but one of the two numerator parallelograms is stretched long and thin, so its area gets small and this would be detected as a false intersection if we didn't check both cases.