/*
* $Id: prune.c,v 1.3 2003/09/11 16:12:48 markus Exp $
*
****************************************************************************
*
* MODULE:       Vector library 
*   	    	
* AUTHOR(S):    Original author CERL, probably Dave Gerdes.
*               Update to GRASS 5.7 Radim Blazek.
*
* PURPOSE:      Lower level functions for reading/writing/manipulating vectors.
*
* COPYRIGHT:    (C) 2001 by the GRASS Development Team
*
*               This program is free software under the GNU General Public
*   	    	License (>=v2). Read the file COPYING that comes with GRASS
*   	    	for details.
*
*****************************************************************************/
/*  @(#)prune.c 3.0  2/19/98  */
/* by Michel Wurtz for GRASS 4.2.1 - <mw@engees.u-strasbg.fr>
 * This is a complete rewriting of the previous dig_prune subroutine.
 * The goal remains : it resamples a dense string of x,y coordinates to
 * produce a set of coordinates that approaches hand digitizing.
 * That is, the density of points is very low on straight lines, and
 * highest on tight curves.
 *
 * The algorithm used is very different, and based on the suppression
 * of intermediate points, when they are closer than thresh from a
 * moving straight line.  
 *
 * The distance between a point M                ->   ->
 * and a AD segment is given                  || AM ^ AD ||
 * by the following formula :            d = ---------------
 *                                                  ->
 *                                               || AD ||
 *
 *                     ->   ->                             ->
 * When comparing   || AM ^ AD ||   and    t = thresh * || AD ||
 *
 *                     ->   ->       ->   ->
 * we call  sqdist = | AM ^ AD | = | OA ^ OM + beta | 
 *
 *                  ->   ->
 *  with     beta = OA ^ OD 
 *
 * The implementation is based on an old integer routine (optimised
 * for machine without math coprocessor), itself inspired by a PL/1
 * routine written after a fortran programm on some prehistoric
 * hardware (IBM 360 probably). Yeah, I'm older than before :-)
 *
 * The algorithm used doesn't eliminate "duplicate" points (following
 * points with same coordinates).  So we should clean the set of points
 * before.  As a side effect, dig_prune can be called with a null thresh
 * value.  In this case only cleaning is made. The command v.prune is to
 * be modified accordingly.
 *
 * Another important note : Don't try too big threshold, this subroutine
 * may produce strange things with some pattern (mainly loops, or crossing
 * of level curves): Try the set { 0,0 -5,0 -4,10 -6,20 -5,30 -5,20 10,10}
 * with a thershold of 5. This isn't a programmation, but a conceptal bug ;-)
 *
 * Input parameters :
 * points->x, ->y   - double precision sets of coordinates.
 * points->n_points - the total number of points in the set.
 * thresh           - the distance that a string must wander from a straight
 *                    line before another point is selected.
 *
 * Value returned : - the final number of points in the pruned set.
 */

#include <stdio.h>
#include "Vect.h"
#include <math.h>

int 
dig_prune (struct line_pnts *points, double thresh)
{
  double *ox, *oy, *nx, *ny;
  double cur_x, cur_y;
  int o_num;
  int n_num;			/* points left */
  int at_num;
  int ij,			/* position of farest point */
    ja, jd, i, j, k, n, inu, it;	/* indicateur de parcours du segment */

  double sqdist;		/* square of distance */
  double fpdist;		/* square of distance from chord to farest point */
  double t, beta;		/* as explained in commented algorithm  */

  double dx, dy;		/* temporary variables */

  double sx[18], sy[18];	/* temporary table for processing points */
  int nt[17], nu[17];

  /* nothing to do if less than 3 points ! */
  if (points->n_points <= 2)
    return (points->n_points);

  ox = points->x;
  oy = points->y;
  nx = points->x;
  ny = points->y;

  o_num = points->n_points;
  n_num = 0;

  /* Eliminate duplicate points */

  at_num = 0;
  while (at_num < o_num)
    {
      if (nx != ox)
	{
	  *nx = *ox++;
	  *ny = *oy++;
	}
      else
	{
	  ox++;
	  oy++;
	}
      cur_x = *nx++;
      cur_y = *ny++;
      n_num++;
      at_num++;

      while (*ox == cur_x && *oy == cur_y)
	{
	  if (at_num == o_num)
	    break;
	  at_num++;
	  ox++;
	  oy++;
	}
    }

  /*  Return if less than 3 points left.  When all points are identical,
   *  output only one point (is this valid for calling function ?) */

  if (n_num <= 2)
    return n_num;

  if (thresh == 0.0)		/* Thresh is null, nothing more to do */
    return n_num;

  /* some (re)initialisations */

  o_num = n_num;
  ox = points->x;
  oy = points->y;

  sx[0] = ox[0];
  sy[0] = oy[0];
  n_num = 1;
  at_num = 2;
  k = 1;
  sx[1] = ox[1];
  sy[1] = oy[1];
  nu[0] = 9;
  nu[1] = 0;
  inu = 2;

  while (at_num < o_num)
    {				/* Position of last point to be    */
      if (o_num - at_num > 14)	/* processed in a segment.         */
	n = at_num + 9;		/* There must be at least 6 points */
      else			/* in the current segment.         */
	n = o_num;
      sx[0] = sx[nu[1]];	/* Last point written becomes      */
      sy[0] = sy[nu[1]];	/* first of new segment.           */
      if (inu > 1)		/* One point was keeped in the     */
	{			/* previous segment :              */
	  sx[1] = sx[k];	/* Last point of the old segment   */
	  sy[1] = sy[k];	/* becomes second of the new one.  */
	  k = 1;
	}
      else
	{			/* No point keeped : farest point  */
	  sx[1] = sx[ij];	/* is loaded in second position    */
	  sy[1] = sy[ij];	/* to avoid cutting lines with     */
	  sx[2] = sx[k];	/* small cuvature.                 */
	  sy[2] = sy[k];	/* First point of previous segment */
	  k = 2;		/* becomes the third one.          */
	}
      /* Loding remaining points         */
      for (j = at_num; j < n; j++)
	{
	  k++;
	  sx[k] = ox[j];
	  sy[k] = oy[j];
	}

      jd = 0;
      ja = k;
      nt[0] = 0;
      nu[0] = k;
      inu = 0;
      it = 0;
      for (;;)
	{
	  if (jd + 1 == ja)	/* Exploration of segment terminated */
	    goto endseg;

	  dx = sx[ja] - sx[jd];
	  dy = sy[ja] - sy[jd];
	  t = thresh * hypot (dx, dy);
	  beta = sx[jd] * sy[ja] - sx[ja] * sy[jd];

	  /* Initializing ij, we don't take 0 as initial value
	     ** for fpdist, in case ja and jd are the same
	   */
	  ij = (ja + jd + 1) >> 1;
	  fpdist = 1.0;

	  for (j = jd + 1; j < ja; j++)
	    {
	      sqdist = fabs (dx * sy[j] - dy * sx[j] + beta);
	      if (sqdist > fpdist)
		{
		  ij = j;
		  fpdist = sqdist;
		}
	    }
	  if (fpdist > t)	/* We found a point to be keeped    */
	    {			/* Restart from farest point        */
	      jd = ij;
	      nt[++it] = ij;
	    }
	  else
	  endseg:{		/* All points are inside threshold. */
	      /* Former start becomes new end     */
	      nu[++inu] = jd;
	      if (--it < 0)
		break;
	      ja = jd;
	      jd = nt[it];
	    }
	}
      for (j = inu - 1; j > 0; j--)
	{			/* Copy of segment's keeped points  */
	  i = nu[j];
	  ox[n_num] = sx[i];
	  oy[n_num] = sy[i];
	  n_num++;
	}
      at_num = n;
    }
  i = nu[0];
  ox[n_num] = sx[i];
  oy[n_num] = sy[i];
  n_num++;
  return n_num;
}