/******************************************************************************
  These are user-defined functions that can be bound to a Postgres backend
  and called by Postgres to execute SQL functions of the same name.

  The calling format for these functions is defined by the CREATE FUNCTION
  SQL statement that binds them to the backend.
********************************************
* Alex Shevlakov, 04/2002 sixote@yahoo.com *
*****************************************************************************/

#include "postgres.h"			/* general Postgres declarations */

#include "fmgr.h"				/* for argument/result macros */
#include "executor/executor.h"	/* for GetAttributeByName() */
#include "utils/geo_decls.h"	/* for point type */
#include <stdio.h>

#ifndef PI
#define PI 3.1415926536
#endif

/* These prototypes just prevent possible warnings from gcc. */
/*check if a point is within buffer of line*/
Datum		path_buffer_contain_pt(PG_FUNCTION_ARGS);
/*removes multiple repeated points from path*/ 
Datum path_without_doubles (PG_FUNCTION_ARGS);
static  PATH *		path_no_dbles(PATH *path);
/*returns closed path which is buffer to another path*/
Datum		return_path_buffer(PG_FUNCTION_ARGS);
/*used to insert points along circle segment between two ends of non-intersecting buffer segments, smoothing the buffer*/
static Point * point_rotate (Point *lpoint, Point *ax_point, float ang,
int k, int n, int napr);
/*writes path to arcinfo UNGENERATE format which is particularly useful for viewing buffer in GRASS*/
Datum write_path_to_file(PG_FUNCTION_ARGS);
int write_path_to_file_internal(PATH *path, char *str );
/*generalize path,i.e., leave each third, or fourth, etc., point*/
Datum reduce_path_points(PG_FUNCTION_ARGS);
static PATH 	*path_reduce(PATH * path, int n_reduce);


/* By Value */


PG_FUNCTION_INFO_V1(path_buffer_contain_pt);

Datum
path_buffer_contain_pt(PG_FUNCTION_ARGS)
{
	PATH	   *path = PG_GETARG_PATH_P(0);
	Point	   *point = PG_GETARG_POINT_P(1);
	float8		radius = PG_GETARG_FLOAT8(2);
	
	CIRCLE		*circle;
	Point	   *pl, *pr, *p_img, *factor;
	double		d;
	int		i;
	

	
	if (radius <= 0) {
		elog(ERROR, "Radius must be a positive number");
		PG_RETURN_NULL();
	}

/* first we check if the circles around vertices contain the point */

	circle = (CIRCLE *) palloc(sizeof(CIRCLE));
	circle->radius = radius;
	
	path =  path_no_dbles(path);
	
	for (i = 0; i < path->npts; i++)
	{
		
		circle->center.x = path->p[i].x;
		circle->center.y = path->p[i].y;

		d = point_dt(&circle->center, point);
		if (d <= circle->radius) {
			pfree(circle);
			PG_RETURN_BOOL(1);
		}		
	}
	
	pfree(circle);

/* now check if the rotated rectangles contain the point */
	
	factor = (Point *) palloc(sizeof(Point));
	
	for (i = 0; i < path->npts-1; i++)
	{
		d = point_dt(&path->p[i+1], &path->p[i]);
		

		
		factor->x = (path->p[i+1].x - path->p[i].x)/d;
		factor->y = (-1) * (path->p[i+1].y - path->p[i].y)/d;
		

		pl = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i]),
						PointPGetDatum(factor)));
		pr = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i+1]),
						PointPGetDatum(factor)));
		p_img = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(point),
						PointPGetDatum(factor)));
		
		if ((p_img->x < pr->x) && (p_img->x > pl->x) 
		&& (p_img->y < (pr->y + radius)) && (p_img->y > (pr->y - radius)) )  {
			
			pfree(pl);
			pfree(pr);
			pfree(p_img);		
			pfree(factor);
			PG_RETURN_BOOL(1);
		}
		
		pfree(pl);
		pfree(pr);
		pfree(p_img);		

		
	}
	
	pfree(factor);

	PG_RETURN_BOOL(0);
}

/*deletes chains of equal points leaving one */

PG_FUNCTION_INFO_V1(path_without_doubles);

Datum 
path_without_doubles (PG_FUNCTION_ARGS)
{
	PATH	   *path = PG_GETARG_PATH_P(0);
	
	PATH *result;
	result = path_no_dbles(path);
	PG_RETURN_PATH_P(result);
}
	
static PATH *
path_no_dbles(PATH * path) {
	
	PATH 		*result;
	int size,i, counter=0;
	Point *ptr;
	
	for (i=0; i < (path->npts) - 1; i++) {
		if ( (path->p[i+1].x == path->p[i].x) &&
		(path->p[i+1].y == path->p[i].y)) counter++;
			
		
	}

	size = offsetof(PATH, p[0]) + sizeof(path->p[0]) * ((path->npts) - counter);

	result = (PATH *) palloc(size);
	
	result->closed = path->closed;
	result->size = size;
	result->npts = (path->npts) - counter;

	ptr  = path->p;
		
	for (i=0; i < result->npts; i++) {
		

		while(1){
		if ( (ptr[0].x == ptr[1].x) &&
		(ptr[0].y == ptr[1].y)) ++ptr;
		else break;
		} 
		
		result->p[i].x = ptr[0].x;
		result->p[i].y = ptr[0].y;
		++ptr;

	}
	
	return result;
	
}

PG_FUNCTION_INFO_V1(return_path_buffer);

Datum
return_path_buffer(PG_FUNCTION_ARGS) {

	PATH	   *path = PG_GETARG_PATH_P(0);
	float8		radius = PG_GETARG_FLOAT8(1);
	
	PATH	   *result; /*for whole buffer line*/
	PATH	   *upper; /*the upper buffer line*/
	PATH	   *lower; /*the lower buffer line*/

	Point	   *pl1, *pr1, *pl2, *pr2; /*images of the two segments*/
	Point	   *pld1, *prd1, *pld2, *prd2; /*down buffer segments*/
	Point	   *plu1, *pru1, *plu2, *pru2; /*upper buffer segments*/
		
	Point	    *factor1, *factor2; /*rotation factor*/
	double		dist1, dist2, deltaf, tang,kappa;
	
	int			size;
	int			i,k, np,n=3;
	int 			lcnt=0, ucnt=0;

		
	/* n - number of vertices in arc 
	connecting two non-intersecting segments of buffer*/

	if((path->p[0].x == path->p[path->npts-1].x) && 
		(path->p[0].y == path->p[path->npts-1].y)) 
		--path->npts;
	

	path =  path_no_dbles(path);
	
	np = path->npts;
	
	factor1 = (Point *) palloc(sizeof(Point));
	factor2 = (Point *) palloc(sizeof(Point));
	plu1 = (Point *) palloc(sizeof(Point));
	pru1 = (Point *) palloc(sizeof(Point));
	pld1 = (Point *) palloc(sizeof(Point));
	prd1 = (Point *) palloc(sizeof(Point));
	plu2 = (Point *) palloc(sizeof(Point));
	pru2 = (Point *) palloc(sizeof(Point));
	pld2 = (Point *) palloc(sizeof(Point));
	prd2 = (Point *) palloc(sizeof(Point));
	
	size = offsetof(PATH, p[0]) +sizeof(path->p[0]) * (np + (np-2)*(n+1) + 1);
	upper = (PATH *) palloc(size);
	lower = (PATH *) palloc(size);
	upper->size = size;
	upper->npts = (np + (np-2)*(n+1)+1);
	upper->closed = path->closed;
	lower->size = size;
	lower->npts = (np + (np-2)*(n+1)+1);
	lower->closed = path->closed;
	
		
   for (i = 0; i < np-2; i++) {
	
	dist1 = point_dt(&path->p[i+1], &path->p[i]);
	
	factor1->x = (path->p[i+1].x - path->p[i].x)/dist1;
	factor1->y = (-1) *(path->p[i+1].y - path->p[i].y)/dist1;

/*first turning them*/	
	pl1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i]),
						PointPGetDatum(factor1)));
	pr1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i+1]),
						PointPGetDatum(factor1)));



/*lifting points above line*/
	plu1->x = pl1->x;
	plu1->y = pl1->y + radius;
	pru1->x = pr1->x;
	pru1->y = pr1->y + radius;
/*lowering points under line*/
	pld1->x = pl1->x;
	pld1->y = pl1->y - radius;
	prd1->x = pr1->x;
	prd1->y = pr1->y - radius;
/*turning them back*/

	factor1->y *= (-1);
	plu1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(plu1),
						PointPGetDatum(factor1)));
	pru1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(pru1),
						PointPGetDatum(factor1)));
	pld1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(pld1),
						PointPGetDatum(factor1)));
	prd1 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(prd1),
						PointPGetDatum(factor1)));


	factor1->y *= (-1);
	

	dist2 = point_dt(&path->p[i+2], &path->p[i+1]);

	factor2->x = (path->p[i+2].x - path->p[i+1].x)/dist2;
	factor2->y = (-1) * (path->p[i+2].y - path->p[i+1].y)/dist2;
	
	pl2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i+1]),
						PointPGetDatum(factor2)));
	pr2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(&path->p[i+2]),
/*lifting points above line*/			PointPGetDatum(factor2)));
	plu2->x = pl2->x;
	plu2->y = pl2->y + radius;
	pru2->x = pr2->x;
	pru2->y = pr2->y + radius;
/*lowering points under line*/
	pld2->x = pl2->x;
	pld2->y = pl2->y - radius;
	prd2->x = pr2->x;
	prd2->y = pr2->y - radius;
/*turning them back*/
	
	factor2->y *= (-1);
	plu2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(plu2),
						PointPGetDatum(factor2)));
	pru2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(pru2),
						PointPGetDatum(factor2)));
	pld2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(pld2),
						PointPGetDatum(factor2)));
	prd2 = DatumGetPointP(DirectFunctionCall2(point_mul,
						PointPGetDatum(prd2),
						PointPGetDatum(factor2)));


	factor2->y *= (-1);

/*lets's compare factors:
sin(b-a) = sin(b)*cos(a) - sin(a)cos(b) >0,<0,==0*/

	deltaf =  (factor1->y)*(factor2->x) - (factor2->y)*(factor1->x);
	
/*tang(phi/2) = mod(deltaf)/1+cos(a-b)*/

	tang= fabs(deltaf)/(1+sqrt(1-deltaf*deltaf));
	kappa = radius * tang/dist1;
	

if (!i) {
	upper->p[i].x=plu1->x;
	upper->p[i].y=plu1->y;
	lower->p[i].x=pld1->x;
	lower->p[i].y=pld1->y;
}

	if (!deltaf) {
	
		upper->p[i+ucnt*(n+1)+1].x=pru1->x;
		upper->p[i+ucnt*(n+1)+1].y=pru1->y;
		upper->p[i+ucnt*(n+1)+2].x=pru2->x;
		upper->p[i+ucnt*(n+1)+2].y=pru2->y;
	
		lower->p[i+1+lcnt*(n+1)].x=prd1->x;
		lower->p[i+1+lcnt*(n+1)].y=prd1->y;
		lower->p[i+2+lcnt*(n+1)].x=prd2->x;
		lower->p[i+2+lcnt*(n+1)].y=prd2->y;
		


	} else if (deltaf > 0) {
/*deltaf >0 - (b-a)>0*/
		
		upper->p[i+ucnt*(n+1)+1].x=pru1->x - (pru1->x - plu1->x) * kappa;
		upper->p[i+ucnt*(n+1)+1].y=pru1->y - (pru1->y - plu1->y) * kappa;
		upper->p[i+ucnt*(n+1)+2].x=pru2->x;
		upper->p[i+ucnt*(n+1)+2].y=pru2->y;
		
		lower->p[i+1+lcnt*(n+1)].x=prd1->x;
		lower->p[i+1+lcnt*(n+1)].y=prd1->y;
		
		
		for (k=1; k<=n; k++) {
		lower->p[i+1+lcnt*(n+1)+k].x =
		(point_rotate(prd1,&path->p[i+1],2*(atan(tang)),k,n,-1))->x;
		lower->p[i+1+lcnt*(n+1)+k].y = 
		(point_rotate(prd1,&path->p[i+1],2*(atan(tang)),k,n,-1))->y;		
		}
		++lcnt;
		
		lower->p[i+1+lcnt*(n+1)].x=pld2->x;
		lower->p[i+1+lcnt*(n+1)].y=pld2->y;
		lower->p[i+2+lcnt*(n+1)].x=prd2->x;
		lower->p[i+2+lcnt*(n+1)].y=prd2->y;
		


	} else {
		
/*deltaf <0 - (b-a)<0*/
		
		upper->p[i+ucnt*(n+1)+1].x=pru1->x;
		upper->p[i+ucnt*(n+1)+1].y=pru1->y;
		
		for (k=1; k<=n; k++) {
		upper->p[i+ucnt*(n+1)+1+k].x = 
		(point_rotate(pru1,&path->p[i+1],2*(atan(tang)),k,n,1))->x;
		upper->p[i+ucnt*(n+1)+1+k].y = 
		(point_rotate(pru1,&path->p[i+1],2*(atan(tang)),k,n,1))->y;
		
		}

		++ucnt;
		
		upper->p[i+ucnt*(n+1)+1].x=plu2->x;
		upper->p[i+ucnt*(n+1)+1].y=plu2->y;
		upper->p[i+ucnt*(n+1)+2].x=pru2->x;
		upper->p[i+ucnt*(n+1)+2].y=pru2->y;
		
		
		lower->p[i+1+lcnt*(n+1)].x=prd1->x - (prd1->x - pld1->x) * kappa;
		lower->p[i+1+lcnt*(n+1)].y=prd1->y - (prd1->y - pld1->y) * kappa;
		lower->p[i+2+lcnt*(n+1)].x=prd2->x;
		lower->p[i+2+lcnt*(n+1)].y=prd2->y;
		
		
	}

	 

	pfree(pl1);
	pfree(pr1);
	pfree(pl2);
	pfree(pr2);

   }

	upper = path_no_dbles(upper);
	lower = path_no_dbles(lower);

/*trimming off the last zero*/

	--lower->npts;
	--upper->npts;
		
	size = offsetof(PATH, p[0]) +sizeof(path->p[0]) * (lower->npts + upper->npts +
	2*n + 1);
	result = (PATH *) palloc(size);

	result->size = size;
	result->npts = (lower->npts + upper->npts + 2*n+1);
	result->closed = lower->closed;
	
	for (i = 0; i < upper->npts; i++)
	{
		result->p[i].x = upper->p[i].x;
		result->p[i].y = upper->p[i].y;
	}
	for (i = 1; i <= n; i++)
	{
		result->p[i + upper->npts - 1].x = 
		(point_rotate(&upper->p[upper->npts - 1],&path->p[np-1],PI,i,n,1))->x;
		result->p[i + upper->npts - 1].y = 
		(point_rotate(&upper->p[upper->npts - 1],&path->p[np-1],PI,i,n,1))->y;
	}
	for (i = 0; i < lower->npts; i++)
	{
		result->p[i + n + upper->npts].x = lower->p[(lower->npts) -1 -i].x;
		result->p[i + n + upper->npts].y = lower->p[(lower->npts) -1 -i].y;
	}
	for (i = 1; i <= n; i++)
	{
		result->p[i + lower->npts + n + upper->npts - 1].x = 
		(point_rotate(&lower->p[0],&path->p[0],PI,i,n,1))->x;
		result->p[i + lower->npts + n + upper->npts - 1].y = 
		(point_rotate(&lower->p[0],&path->p[0],PI,i,n,1))->y;
	}
	
	result->p[lower->npts + upper->npts + 2*n].x = result->p[0].x;
	result->p[lower->npts + upper->npts + 2*n].y = result->p[0].y;
	


pfree(upper);
pfree(lower);
pfree(factor1);
pfree(factor2);
pfree(plu1);
pfree(pru1);
pfree(pld1);
pfree(prd1);
pfree(plu2);
pfree(pru2);
pfree(pld2);
pfree(prd2);	

PG_RETURN_PATH_P(result);	
}

static Point *
point_rotate (Point *lpoint, Point *ax_point, float ang, int k, int n, int napr) {

Point *factor, *result;

factor = (Point *) palloc(sizeof(Point));
result = (Point *) palloc(sizeof(Point));

if(!napr) 
elog(ERROR,"Parameter of rotation direction can't be zero, 
	must be integer positive or negative)");
napr = abs(napr)/napr;

factor->x = cos((ang/(n+1)) * k);
factor->y = (-1) * sin((ang/(n+1)) * k) * napr;

lpoint = DatumGetPointP(DirectFunctionCall2(point_sub,
					PointPGetDatum(lpoint),
					PointPGetDatum(ax_point)));
result = DatumGetPointP(DirectFunctionCall2(point_mul,
					PointPGetDatum(lpoint),
					PointPGetDatum(factor)));
result = DatumGetPointP(DirectFunctionCall2(point_add,
					PointPGetDatum(result),
					PointPGetDatum(ax_point)));

pfree(factor);
return result;
}

PG_FUNCTION_INFO_V1(write_path_to_file);

Datum 
write_path_to_file(PG_FUNCTION_ARGS) {
PATH	   *path = PG_GETARG_PATH_P(0);

int result;
char *str = "/tmp/boundary.pol";

result = write_path_to_file_internal(path, str );
PG_RETURN_BOOL(result == 0); 
};

int
write_path_to_file_internal(PATH *path, char *str ) {

	
	FILE *fp;
	int i;
	
	fp = fopen(str,"w");
	if ( fp == NULL) return (-1);
	fprintf(fp,"1\n");
	for (i=0;i<path->npts;i++) {
	
	fprintf(fp,"%f %f\n", path->p[i].x, path->p[i].y);
	}
	fprintf(fp,"END\n");
	fprintf(fp,"END\n");
	fclose(fp);
	return 0;
}

PG_FUNCTION_INFO_V1(reduce_path_points);
Datum 
reduce_path_points(PG_FUNCTION_ARGS) {
PATH	   *path = PG_GETARG_PATH_P(0);
int 		n_reduce = PG_GETARG_INT32(1);
PATH *result;
result = path_reduce(path, n_reduce);
PG_RETURN_PATH_P(result);
}

static PATH *
path_reduce(PATH * path, int n_reduce) {
	
	PATH 		*result;
	int size,i,k;
	Point *ptr;


	size = offsetof(PATH, p[0]) + sizeof(path->p[0]) * (floor((path->npts)/n_reduce) + 1);

	result = (PATH *) palloc(size);
	
	result->closed = path->closed;
	result->size = size;
	result->npts = ceil((path->npts)/n_reduce);

	ptr  = path->p;
		
	for (i=0; i < result->npts-1; i++) {
		

		
		result->p[i].x = ptr[0].x;
		result->p[i].y = ptr[0].y;
		for (k=1; k<=n_reduce; k++)
				++ptr;

	}
	result->p[result->npts-1].x = path->p[path->npts-1].x;
	result->p[result->npts-1].y = path->p[path->npts-1].y;
	return result;
	
}