This function is fine for basic polygons but doesn't scale well since its checking all coordinates for every y-pixel.

Heres an optimized version. Basic logic remains the same this just maintains an ordered list of intersections,
tracking in-out points, to avoid re-computing every row, this means sorting is only done when out of order segments are found and once sorted, the segments only need to be re-ordered when the cross each other.

It also makes a minor change, rounding the value when converting to an int,
so the result doesn't change for positive/negative polygons.

Using the lasso tool is a quick way to see its working.
I have some tests here but not in C/C++, ideally this patch would include gtests too.

Simple timing test shows 7.5x speedup (11x if the call to round is removed) for 10k points over a 2560x1580 region.

Also, the function is getting a bit involved and IMHO is outside the scope of BLI_math,
would suggest to move this and to a different module (BLI_raster ?), along with plot_line_v2v2i.

Note, one weakness in this method is that newly intersected nodes are added to the end of the array and need to be re-sorted.
For more efficient handling of dense polygons (100 - 1000+ intersections at once on the Y axis), it might be better to use a b-tree however that would make the patch a lot more involved.
Nevertheless, it improves over sorting every time.