How to solve the lack of symmetry in this line rasterization algorithm?


I am working on an algorithm to snap 3D lines to a 3D tessellated space.

Here is a 2D example of the algorithm that works for positive and negative slopes between -1 and 1 inclusive.

Using the slope to calculate the value of y at each x results in slow and error-prone floating calculations. The solution is to simulate the division with a remainder variable.

When dx >= dy, start with an initial remainder variable ry = 0. Then, for each x increment, add dy to ry variable. When it surpasses dx, increment y, then set ry equal to ry - dx.

function line(x1, y1, x2, y2) {
  let points = []
  let dx = Math.abs(x2 - x1);
  let dy = Math.abs(y2 - y1);
  // The remainder variable for y axes. 
  // No rx is created because we are assuming dx >= dy.
  let ry = 0;
  // Current value of y for a given point
  let y = 0;

  // The slope could be positive or negative, so increments coordinates as they go down or up.
  let pointIncrement;
  if (x2 > x1) {
    pointIncrement = 1;
  } else if (x2 < x1) {
    pointIncrement = -1;
    y = y1
  for (let x = x1; pointIncrement < 0 ? x >= x2 : x <= x2; x += pointIncrement) {
    if (ry >= dx) {
      ry -= dx;
      y += pointIncrement;
    // Add dy to ry until it surpasses dx. This simulates the division of dy/dx for slope.
    ry += dy;
    points.push([x, y])
  return points

Now, if you call the function with a slope of 1/4th:


You get the following results:


Now, if you call it again but in the negative direction, then reverse the coordinate order:


You get the following results:


Why is this occurring?

Is anyone aware of a solution to this problem to make the negative slope symmetric to the positive slope?