/*
 *  avhrr.c -- Operations on AVHRR tuples.
 */

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

#include "tmp/c.h"
#include "tmp/postgres.h"
#include "tmp/libpq-fe.h"

#include "utils/rel.h"
#include "utils/log.h"

#include "access/htup.h"

#include "catalog/syscache.h"
#include "catalog/pg_proc.h"
#include "catalog/pg_type.h"

#include "geoloc.h"
#include "array.h"

/* NEVER define this if you're using the postmaster! */
/* #define	AVHRR_DEBUG */

ObjectId	avhrr_findfunc();
double		geo_max();

/*
 *  clip_avhrr -- Clip an AVHRR image to a bounding rectangle.
 *
 *	This routine converts the input geolocs into LAZEA geolocs, finds
 *	the minimum bounding LAZEA rectangle for the input geoloc rectangle,
 *	and clips the AVHRR image to that rectangle.  A new large object
 *	is created, and an array item is created, filled in, and returned.
 *
 *	The original goal for the demo was to have arrays of different
 *	formats, with functions dispatching through a FORMATS table in
 *	the same way that GEOLOC does.  The deadline was too close; we
 *	hardwire a few formats in the code below, and ask the user to
 *	pass in the target format to the clip routine.
 *
 *	The fmt argument is actually a char8.
 */

ARRAY *
clip_avhrr(avhrr_tup, fmt, nw, se)
	TUPLE *avhrr_tup;
	char *fmt;
	GEOLOC *nw;
	GEOLOC *se;
{
	char *retval;
	GEOLOC *imgnw, *imgse;
	GEOLOC *sw, *ne;
	ARRAY *imgarray;
	ARRAY *result;
	ObjectId fromstd, tostd, clip, outfunc;
	char *s;
	double minx, maxx, miny, maxy;
	int upperleft_x, upperleft_y, lowerright_x, lowerright_y;
	Boolean isnull;

	imgnw = (GEOLOC *) GetAttributeByName(avhrr_tup, "image_nw", &isnull);
	if (isnull)
		elog(WARN, "AVHRR tuple has a NULL image_nw coordinate");

	imgse = (GEOLOC *) GetAttributeByName(avhrr_tup, "image_se", &isnull);
	if (isnull)
		elog(WARN, "AVHRR tuple has a NULL image_se coordinate");

	imgarray = (ARRAY *) GetAttributeByName(avhrr_tup, "image", &isnull);
	if (isnull)
		elog(WARN, "AVHRR tuple has no image large object");

	/* this'll never happen... */
	if (strncmp(&(nw->geo_name[0]), &(se->geo_name[0]), 8) != 0)
		elog(WARN, "clip coordinates must be in the same projection");

	/*
	 *  Construct the southwest and northeast corners of the clip
	 *  rectangle supplied by the user.
	 */

	sw = (GEOLOC *) palloc(nw->geo_varlen);
	bcopy((char *) nw, (char *) sw, nw->geo_varlen);
	sw->geo_y = se->geo_y;

	ne = (GEOLOC *) palloc(se->geo_varlen);
	bcopy((char *) se, (char *) ne, se->geo_varlen);
	ne->geo_y = nw->geo_y;

#ifdef	AVHRR_DEBUG
	outfunc = avhrr_findfunc("geoloc_out", 1);
	s = (char *) fmgr(outfunc, sw);
	fprintf(stderr, "sw\t(%s)\n", s);
	s = (char *) fmgr(outfunc, nw);
	fprintf(stderr, "nw\t(%s)\n", s);
	s = (char *) fmgr(outfunc, ne);
	fprintf(stderr, "ne\t(%s)\n", s);
	s = (char *) fmgr(outfunc, se);
	fprintf(stderr, "se\t(%s)\n", s);
#endif	/* AVHRR_DEBUG */

	/*
	 *  Find the function IDs of the functions we need to do coordinate
	 *  conversions.
	 */

	fromstd = avhrr_findfunc("geoloc_fromstd", 2);
	tostd = avhrr_findfunc("geoloc_tostd", 1);

	/*
	 *  Convert each of the corners of the bounding rectangle from
	 *  the user to LAZEA.  The fmgr invocations below are equivalent
	 *  to the postquel
	 *
	 *	coord = geoloc_fromstd(geoloc_tostd(coord), "lazea")
	 */

	nw = (GEOLOC *) fmgr(fromstd, (GEOLOC *) fmgr(tostd, nw), "lazea");
	sw = (GEOLOC *) fmgr(fromstd, (GEOLOC *) fmgr(tostd, sw), "lazea");
	ne = (GEOLOC *) fmgr(fromstd, (GEOLOC *) fmgr(tostd, ne), "lazea");
	se = (GEOLOC *) fmgr(fromstd, (GEOLOC *) fmgr(tostd, se), "lazea");

#ifdef	AVHRR_DEBUG
	s = (char *) fmgr(outfunc, sw);
	fprintf(stderr, "\nsw\t(%s)\n", s);
	s = (char *) fmgr(outfunc, nw);
	fprintf(stderr, "nw\t(%s)\n", s);
	s = (char *) fmgr(outfunc, ne);
	fprintf(stderr, "ne\t(%s)\n", s);
	s = (char *) fmgr(outfunc, se);
	fprintf(stderr, "se\t(%s)\n", s);
#endif	/* AVHRR_DEBUG */

	/*
	 *  Now we have to compute the minimal bounding box for the
	 *  coordinates in the new projection.  This operation actually
	 *  depends on the projection, but we don't have operators for
	 *  the projection, yet, so we need to hardwire behavior for the
	 *  demo.  The solution is to assume that the larger value in
	 *  magnitude is the larger value, regardless of sign.  This
	 *  won't work if the line on any side of the converted rectangle
	 *  (actually four-sided polygon) spans a zero in the projection,
	 *  like the date line or the equator, but the demo is on the
	 *  conterminous US, so this is a safe way to cheat.
	 */

#if 0
	nw->geo_x = geo_max(nw->geo_x, sw->geo_x);
	nw->geo_y = geo_max(nw->geo_y, ne->geo_y);
	se->geo_x = geo_max(se->geo_x, ne->geo_x);
	se->geo_y = geo_max(se->geo_y, sw->geo_y);
#endif

	minx = nw->geo_x;
	if (sw->geo_x < minx)
		minx = sw->geo_x;
	if (se->geo_x < minx)
		minx = se->geo_x;
	if (ne->geo_x < minx)
		minx = ne->geo_x;

	maxx = nw->geo_x;
	if (sw->geo_x > maxx)
		maxx = sw->geo_x;
	if (se->geo_x > maxx)
		maxx = se->geo_x;
	if (ne->geo_x > maxx)
		maxx = ne->geo_x;

	miny = nw->geo_y;
	if (sw->geo_y < miny)
		miny = sw->geo_y;
	if (se->geo_y < miny)
		miny = se->geo_y;
	if (ne->geo_y < miny)
		miny = ne->geo_y;

	maxy = nw->geo_y;
	if (sw->geo_y > maxy)
		maxy = sw->geo_y;
	if (se->geo_y > maxy)
		maxy = se->geo_y;
	if (ne->geo_y > maxy)
		maxy = ne->geo_y;

	if (imgnw->geo_x < imgse->geo_x) {
		if (minx < imgnw->geo_x)
			minx = imgnw->geo_x;
		if (maxx > imgse->geo_x)
			maxx = imgse->geo_x;
		nw->geo_x = minx;
		se->geo_x = maxx;
	} else {
		if (minx > imgnw->geo_x)
			minx = imgnw->geo_x;
		if (maxx < imgse->geo_x)
			maxx = imgse->geo_x;
		nw->geo_x = maxx;
		se->geo_x = minx;
	}

	if (imgnw->geo_y < imgse->geo_y) {
		if (miny < imgnw->geo_y)
			miny = imgnw->geo_y;
		if (maxy > imgse->geo_y)
			maxy = imgse->geo_y;
		nw->geo_y = miny;
		se->geo_y = maxy;
	} else {
		if (miny > imgnw->geo_y)
			miny = imgnw->geo_y;
		if (maxy < imgse->geo_y)
			maxy = imgse->geo_y;
		nw->geo_y = maxy;
		se->geo_y = miny;
	}

#ifdef	AVHRR_DEBUG
	s = (char *) fmgr(outfunc, nw);
	fprintf(stderr, "\nnw\t(%s)\n", s);
	s = (char *) fmgr(outfunc, se);
	fprintf(stderr, "se\t(%s)\n", s);
#endif	/* AVHRR_DEBUG */

	/* done with these now */
	pfree (ne);
	pfree (sw);

	/*
	 *  Coordinates for the array are relative to the image's origin,
	 *  so convert to that coordinate system.  Also, we know that
	 *  AVHRR pixels are 1km on a side, so we scale down by 1000.
	 */

	upperleft_x = floor((nw->geo_x - imgnw->geo_x) / 1000.0);
	upperleft_y = floor((imgnw->geo_y - nw->geo_y) / 1000.0);
	lowerright_x = ceil((se->geo_x - imgnw->geo_x) / 1000.0);
	lowerright_y = ceil((imgnw->geo_y - se->geo_y) / 1000.0);

#ifdef AVHRR_DEBUG
	elog(NOTICE, "nw = (%lg, %lg)", nw->geo_x, nw->geo_y);
	elog(NOTICE, "se = (%lg, %lg)", se->geo_x, se->geo_y);
#endif

	/* finally, arrays in postgres are zero-based */
	if (upperleft_x > 0)
		upperleft_x--;
	if (upperleft_y > 0)
		upperleft_y--;
	if (lowerright_x > 0)
		lowerright_x--;
	if (lowerright_y > 0)
		lowerright_y--;

#ifdef AVHRR_DEBUG
	elog(NOTICE, "ul = (%d, %d)", upperleft_x, upperleft_y);
	elog(NOTICE, "lr = (%d, %d)", lowerright_x, lowerright_y);
#endif

#ifdef	AVHRR_DEBUG
	fprintf(stderr, "\n(%d, %d) to (%d, %d)\n",
		upperleft_x, upperleft_y, lowerright_x, lowerright_y);
#endif	/* AVHRR_DEBUG */

	clip = avhrr_findfunc("arrayclip", 6);
	result = (ARRAY *) fmgr(clip, imgarray, fmt,
				upperleft_x, upperleft_y,
				lowerright_x, lowerright_y);

	return (result);
}

/*
 *  geo_max() -- Compute the maximal of two geographic half-coord values.
 *
 *	We use this routine to compute the minimal bounding rectangle for
 *	a four-sided polygon after we transform a rectangle into a new
 *	coordinate system.  This operation is actually projection-
 *	dependent, but for the demo (which operates on the conterminous
 *	US, and for which we don't cross any interesting zeroes) we
 *	can cheat in the following way:  If the sign of the two components
 *	is the same, then take the one with the larger magnitude.
 *
 *	We truncate to integer doubles, while we're at it.
 */

double
geo_max(a, b)
	double a;
	double b;
{
	if (a < 0.0 && b < 0.0)
		return (floor(a < b ? a : b));
	else if (a > 0.0 && b > 0.0)
		return (ceil(a > b ? a : b));

	elog(WARN, "projection conversion error: can't compare %lg and %lg",
		a, b);
}

/*
 *  avhrr_findfunc() -- Find a function by function name.
 *
 *	This routine looks up the named function in the pg_proc system
 *	cache and returns its function id, for use by the function manager.
 *	We could supply the actual argument types to the function, but
 *	that's a fair amount of trouble.  Instead, we allow the system
 *	to find the function without any type information.  This only
 *	works if there's exactly one such named function in the catalogs.
 */

ObjectId
avhrr_findfunc(funcname, nargs)
	char *funcname;
	int nargs;
{
	int i;
	ObjectId fid;
	bool retset;
	ObjectId rettype;
	ObjectId *actualtypes;
	ObjectId argtypes[8];

	for (i = 0; i < 8; i++)
		argtypes[i] = UNKNOWNOID;

	if (!func_get_detail(funcname, nargs, argtypes, &fid, &rettype,
			     &retset, &actualtypes))
		elog(WARN, "cannot find function %s in %s", funcname,
			   Name_pg_proc);

	return (fid);
}