/* line.c - Basic line drawing.
 * (C) Copyright 1999 John Halleck
 * All rights reserved
 */

/* Version of August 13th, 1999 */

#include "errors.h"
 /* we use package standard error codes. */

#include "line.h"

/*
 * Draw a line between (x1, y1) and (x2, y2)
 *
 *    The equation of the line is..
 *        |  X  Y  1 |
 *        | x1 y1  1 | = 0
 *        | x2 y2  1 |
 * =>
 *        X*(y1-y2)  -  Y*(x1-x2)  +  x1*y2  -  x2*y1 = 0
 *
 *    If you plug other points into this equation, you get a number other
 *    than zero.  This number is an (unnormalized) distance function.
 *    This means that is distance multiplied by some constant.
 *        (The constant is sqrt((y1-y2)**2+(x1-x2)**2) )
 *    since we are only going to compare magnitudes, we don't care
 *    about the constant, and we don't bother to compute it.
 *
 *    So... The distance (error) function we will use is this (unnormalized)
 *    distance from the line.
 *
 *          distance = X*(y1-y2)  -  Y*(x1-x2)  +  x1*y2  -  x2*y1
 *
 *    This distance function has the nice property that it is positive
 *    on one side of the line, and negative on the other.
 *
 *    As the distance is zero at both endpoints, the cross product
 *    terms are zero there.   As the change as we move from grid point
 *    to grid point can be computed from DX and DY, we do not need
 *    to compute the cross product terms at all. 
 *    So, our simplified initial distance function is:
 * 
 *          distance = X * (y1-y2)  -  Y * (x1-x2)
 * 
 * The idea of the line drawing algorithm is that for each step of the line,
 * one wants to pick the direction that takes one closest to the line being
 * drawn.  This is equivalent, in the end, in trying to push us always
 * towards the line.  Since the distance function is signed, the rule
 * reduces to going in the direction that makes it more positive if it
 * is negative, and in the direction that makes it more negative if it
 * is positive.   This is a very clean test in most languages, and the
 * underlying hardware does it well on most machines.
 *
 * If it is the case that both directions would take you the same
 * distance from the line,  you want to pick one in such a manner that
 * the same one is taken for each direction the line could be drawn.
 *
 * If you examine the algorithm carefully, you will note that the
 * choices of direction are either along an axis, or along a diagonal.
 * Moving in the one direction makes the distance more positive, and
 * moving in the other direction makes the distance more negative.
 * Based on the sign of the current distance, we always move in the
 * direction that makes the sign go more the other direction.
 *
 * There is also a version of this algorithm (that produces the same
 * results) based on looking at the result with modular arithmetic.
 * [Derived by LeRoy Eide <LeRoy.Eide@cc.utah.edu>] It can be used
 * to keep a consistent dot pattern across strips when you have to
 * make a plot in sections.   It has a somewhat different setup, and
 * the same loop.
 */ 

/* External support routine needed */
extern error drawpoint (int x, int y);

error line (int x1, int y1,  int x2, int y2)  {
  int axis, axisX, axisY, diag, diagX, diagY, bias;
  /* How far we step, and in what direction, for a step */
  int dist, cX, cY, dX, dY;
  /* Unnormalized distance */
  int numpnt, X, Y;
  /* Current location */
  int Errors;
  /* Any problems we've had. */

  /* Set up to know what direction we are going. */
  dX = x2 - x1; if (dX<0) { cX = -dX; diagX = -1; }
                else      { cX =  dX; diagX =  1; }
  dY = y2 - y1; if (dY<0) { cY = -dY; diagY = -1; }
                else      { cY =  dY; diagY =  1; }

  /* Decide which axis we are using, make the other the diagonal */
  if (cY <= cX) { axisX = diagX; axisY = 0;     axis = cY; bias = cX; }
  else          { axisX = 0;     axisY = diagY; axis = cX; bias = cY; }

  /* Compute up the amount of change going each direction makes. */
  diag = axis - bias; dist = - bias - 1;

/* This line, see: Tran Thong, Tektronix Laboratories, September 19, 1980
 * This line preserves the Tran Thong property in the resultant dots.
 * [Makes it chose exactly the same dots for lines drawn backwards.]
 * (Choice of diagX or diagY was arbitrary)
 * This line can be omitted without any change in line quality,
 * although I can't imagine why you'd want to get rid of it.
 */ if (diagX < 0) dist --;

    /* Bias line so that discontinuities are towards the middle 
     * Imagine a line 20 in x and 1 in y.  Somewhere it has to
     * make the leap between the one y value and the other.
     * This has to be near the middle if the correct properties
     * are to be observed.
     */
   dist = axis + (dist >> 1);

   /* We have to start somewhere. */
   numpnt = bias + 1; X = x1; Y = y1;

   /* ****************************************** */

   /* FINALLY, here is the core loop. */
   while (numpnt-- > 0) {

      /* Put out the current point */
      if ((Errors = drawpoint (X, Y))) return Errors;

      /* Decide where to move, and update the distance function,
       * and the coordinates.
       */
      if (dist < 0 ) { dist += axis; X += axisX; Y += axisY; }
      else           { dist += diag; X += diagX; Y += diagY; }
   }
   return NoError;
}

/* -------------------------------------------------------------------------- */
/* ========================================================================== */
/* -------------------------------------------------------------------------- */

/* The observation has also been made that the dist (error) function
 * above can be used to make anti-aliased lines directly, by distributing
 * the line across the pixel's that straddle it.
 *
 * Basicly, if the line's error value is as far one way as it can be, 
 * then the pixel on that side of the line is as far from the line
 * as it can get.  It should get none of the line, and the other side
 * should get all.   The same argument works the other way.
 * If it is, for example, on the line, the pixels on both sides should
 * get the line's value equally.  And so on.
 *
 * That leads to the code below:
 */

/* Support routine needed. */
extern error drawspoint (int x, int y, double value);

error sline (int x1, int y1, int x2, int y2) {
/* (Comments and structure as above...) */
  int axis, axisX, axisY, diag, diagX, diagY, bias;
  int dist, cX, cY, dX, dY;
  int numpnt, X, Y;
  error Errors;
  double shade, scale;

  /* This setup is the same as the previous routine, I'm not going to
   * repeat all that...  But I'll comment places where they are different.
   */
  dX = x2 - x1; if (dX<0) { cX = -dX; diagX = -1; }
                else      { cX =  dX; diagX =  1; }
  dY = y2 - y1; if (dY<0) { cY = -dY; diagY = -1; }
                else      { cY =  dY; diagY =  1; }

  if (cY <= cX) { axisX = diagX; axisY = 0;     axis = cY; bias = cX; }
  else          { axisX = 0;     axisY = diagY; axis = cX; bias = cY; }

  diag = axis - bias; dist = - bias - 1;

   if (diagX < 0) dist --;

   /* The following is not needed (and in fact wrong) for anti-aliased
    * lines.  So it is commented out here.
    */
   /* dist = axis + (dist >> 1); */ /* anti aliased lines have no jumps */

   /* The scale value is the inverse of the magnitude of the possible range
    * of values that the distance can take.
    */
   scale = 1.0 / (float) (bias+1+axis);
   numpnt = bias + 1; X = x1; Y = y1;

   /* Core loop. */
   while (numpnt-- > 0) {
      shade = scale * (float) (dist+bias+1);
      if (shade != 1.0) {
          Errors = drawspoint (X,       Y,       1.0 - shade);
          if (Errors) return Errors;
      }
      if (shade != 0.0) {
          Errors = drawspoint (X+axisX, Y+axisY,       shade);
          if (Errors) return Errors;
      }
      if (dist < 0 ) { dist += axis; X += axisX; Y += axisY; }
      else           { dist += diag; X += diagX; Y += diagY; }
   }
   return NoError;
}