Re: Line intersection point is wrong

Emre Hasegeli <emre@hasegeli.com> writes:
>> You haven't done anything to exclude the possibility that l1->B is zero,
>> so you could be getting zero-divide in the y calculation.

> Ah, yes.  We can use l2 when l1->B is zero.

After working out the algebra by hand, I think the attached is correct
(and it does pass the regression tests, yay).  I also made the
line_parallel and line_perp tests more symmetric and zero-divide-free.

            regards, tom lane

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index 657bcee..ce796c7 100644
*** a/src/backend/utils/adt/geo_ops.c
--- b/src/backend/utils/adt/geo_ops.c
*************** line_parallel(PG_FUNCTION_ARGS)
*** 1108,1117 ****
      LINE       *l1 = PG_GETARG_LINE_P(0);
      LINE       *l2 = PG_GETARG_LINE_P(1);

!     if (FPzero(l1->B))
          PG_RETURN_BOOL(FPzero(l2->B));

!     PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
  }

  Datum
--- 1108,1119 ----
      LINE       *l1 = PG_GETARG_LINE_P(0);
      LINE       *l2 = PG_GETARG_LINE_P(1);

!     if (FPzero(l1->A))            /* horizontal? */
!         PG_RETURN_BOOL(FPzero(l2->A));
!     if (FPzero(l1->B))            /* vertical? */
          PG_RETURN_BOOL(FPzero(l2->B));

!     PG_RETURN_BOOL(FPeq(l2->A * l1->B, l1->A * l2->B));
  }

  Datum
*************** line_perp(PG_FUNCTION_ARGS)
*** 1120,1131 ****
      LINE       *l1 = PG_GETARG_LINE_P(0);
      LINE       *l2 = PG_GETARG_LINE_P(1);

!     if (FPzero(l1->A))
          PG_RETURN_BOOL(FPzero(l2->B));
!     else if (FPzero(l1->B))
          PG_RETURN_BOOL(FPzero(l2->A));

!     PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
  }

  Datum
--- 1122,1133 ----
      LINE       *l1 = PG_GETARG_LINE_P(0);
      LINE       *l2 = PG_GETARG_LINE_P(1);

!     if (FPzero(l1->A))            /* horizontal? */
          PG_RETURN_BOOL(FPzero(l2->B));
!     if (FPzero(l1->B))            /* vertical? */
          PG_RETURN_BOOL(FPzero(l2->A));

!     PG_RETURN_BOOL(FPeq(l2->A * l1->B, -(l1->A * l2->B)));
  }

  Datum
*************** line_interpt_internal(LINE *l1, LINE *l2
*** 1233,1250 ****

      if (FPzero(l1->B))            /* l1 vertical? */
      {
!         x = l1->C;
!         y = (l2->A * x + l2->C);
      }
      else if (FPzero(l2->B))        /* l2 vertical? */
      {
!         x = l2->C;
!         y = (l1->A * x + l1->C);
      }
      else
      {
!         x = (l1->C - l2->C) / (l2->A - l1->A);
!         y = (l1->A * x + l1->C);
      }
      result = point_construct(x, y);

--- 1235,1252 ----

      if (FPzero(l1->B))            /* l1 vertical? */
      {
!         x = -(l1->C / l1->A);
!         y = -(l2->A * x + l2->C) / l2->B;
      }
      else if (FPzero(l2->B))        /* l2 vertical? */
      {
!         x = -(l2->C / l2->A);
!         y = -(l1->A * x + l1->C) / l1->B;
      }
      else
      {
!         x = (l1->B * l2->C - l2->B * l1->C) / (l1->A * l2->B - l2->A * l1->B);
!         y = -(l1->A * x + l1->C) / l1->B;
      }
      result = point_construct(x, y);