/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * Copyright (C) 2001-2003 Refractions Research Inc.
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU General Public Licence. See the COPYING file.
 *
 **********************************************************************/


/* PostgreSQL headers */
#include "postgres.h"
#include "fmgr.h"
#include "miscadmin.h"
#include "utils/memutils.h"
#include "executor/spi.h"
#include "access/hash.h"
#include "utils/hsearch.h"

/* PostGIS headers */
#include "../postgis_config.h"
#include "liblwgeom.h"
#include "lwgeom_pg.h"
#include "lwgeom_cache.h"
#include "lwgeom_transform.h"

/* C headers */
#include <float.h>
#include <string.h>
#include <stdio.h>

#define maxprojlen  512
#define spibufferlen 512


/*
 * PROJ 4 backend hash table initial hash size
 * (since 16 is the default portal hash table size, and we would
 * typically have 2 entries per portal
 * then we shall use a default size of 256)
 */
#define PROJ_BACKEND_HASH_SIZE 256

/* Global to hold the Proj object cache */
PROJSRSCache *PROJ_CACHE = NULL;


/**
 * Utility structure to get many potential string representations
 * from spatial_ref_sys query.
 */
typedef struct {
	char* authtext; /* auth_name:auth_srid */
	char* srtext;
	char* proj4text;
} PjStrs;

/* Internal Cache API */
static LWPROJ *
AddToPROJSRSCache(PROJSRSCache *PROJCache, int32_t srid_from, int32_t srid_to);
static void DeleteFromPROJSRSCache(PROJSRSCache *PROJCache, uint32_t position);

static void
PROJSRSDestroyPJ(void *projection)
{
	LWPROJ *pj = (LWPROJ *)projection;

	if (pj->pj)
	{
		proj_destroy(pj->pj);
		pj->pj = NULL;
	}
}

/**
* Destroy all the malloc'ed proj objects stored
* in the PROJSRSCache, in case the containing
* context is removed (shouldn't actually happen)
*/
static void
PROJSRSDestroyPortalCache(void *portalCache)
{
	PROJSRSCache* cache = (PROJSRSCache*)portalCache;
	for (uint32_t i = 0; i < cache->PROJSRSCacheCount; i++)
	{
		if (cache->PROJSRSCache[i].projection)
			PROJSRSDestroyPJ(cache->PROJSRSCache[i].projection);
	}
}

/**
* Get the Proj cache entry from the global variable if one exists.
* If it doesn't exist, make a new blank one and return it.
*/
PROJSRSCache *
GetPROJSRSCache()
{
	PROJSRSCache* cache = PROJ_CACHE;
	if (!cache)
	{
		/* Put proj cache in a child of the CacheContext */
		MemoryContext context = AllocSetContextCreate(
		    CacheMemoryContext,
		    "Proj Context",
		    ALLOCSET_SMALL_SIZES);

		/* Allocate in the upper context */
		cache = MemoryContextAllocZero(context, sizeof(PROJSRSCache));

		if (!cache)
			elog(ERROR, "Unable to allocate space for PROJSRSCache in context %p", (void *)context);

		cache->PROJSRSCacheCount = 0;
		cache->PROJSRSCacheContext = context;

		/* Use this to clean up PROJSRSCache in event of MemoryContext reset */
		MemoryContextCallback* callback = MemoryContextAlloc(context, sizeof(MemoryContextCallback));
		callback->func = PROJSRSDestroyPortalCache;
		callback->arg = (void *)cache;
		MemoryContextRegisterResetCallback(context, callback);

		PROJ_CACHE = cache;
	}
	return cache;
}


/*****************************************************************************
 * Per-cache management functions
 */

static LWPROJ *
GetProjectionFromPROJCache(PROJSRSCache *cache, int32_t srid_from, int32_t srid_to)
{
	uint32_t i;
	for (i = 0; i < cache->PROJSRSCacheCount; i++)
	{
		if (cache->PROJSRSCache[i].srid_from == srid_from &&
		    cache->PROJSRSCache[i].srid_to == srid_to)
		{
			cache->PROJSRSCache[i].hits++;
			return cache->PROJSRSCache[i].projection;
		}
	}

	return NULL;
}

static char*
SPI_pstrdup(const char* str)
{
	char* ostr = NULL;
	if (str)
	{
		ostr = SPI_palloc(strlen(str)+1);
		strcpy(ostr, str);
	}
	return ostr;
}

static PjStrs
GetProjStringsSPI(int32_t srid)
{
	int spi_result;
	char proj_spi_buffer[spibufferlen];
	PjStrs strs;
	memset(&strs, 0, sizeof(strs));

	/* Connect */
	spi_result = SPI_connect();
	if (spi_result != SPI_OK_CONNECT)
	{
		elog(ERROR, "Could not connect to database using SPI");
	}

	static char *proj_str_tmpl =
	    "SELECT proj4text, auth_name, auth_srid, srtext "
	    "FROM %s "
	    "WHERE srid = %d "
	    "LIMIT 1";
	snprintf(proj_spi_buffer, spibufferlen, proj_str_tmpl, postgis_spatial_ref_sys(), srid);

	/* Execute the query, noting the readonly status of this SQL */
	spi_result = SPI_execute(proj_spi_buffer, true, 1);

	/* Read back the PROJ text */
	if (spi_result == SPI_OK_SELECT && SPI_processed > 0)
	{
		/* Select the first (and only tuple) */
		TupleDesc tupdesc = SPI_tuptable->tupdesc;
		SPITupleTable *tuptable = SPI_tuptable;
		HeapTuple tuple = tuptable->vals[0];
		/* Always return the proj4text */
		char* proj4text = SPI_getvalue(tuple, tupdesc, 1);
		if (proj4text && strlen(proj4text))
			strs.proj4text = SPI_pstrdup(proj4text);

		/* For Proj >= 6 prefer "AUTHNAME:AUTHSRID" to proj strings */
		/* as proj_create_crs_to_crs() will give us more consistent */
		/* results with authority numbers than with proj strings */
		char* authname = SPI_getvalue(tuple, tupdesc, 2);
		char* authsrid = SPI_getvalue(tuple, tupdesc, 3);
		if (authname && authsrid &&
		    strlen(authname) && strlen(authsrid))
		{
			char tmp[maxprojlen];
			snprintf(tmp, maxprojlen, "%s:%s", authname, authsrid);
			strs.authtext = SPI_pstrdup(tmp);
		}

		/* Proj6+ can parse srtext, so return that too */
		char* srtext = SPI_getvalue(tuple, tupdesc, 4);
		if (srtext && strlen(srtext))
			strs.srtext = SPI_pstrdup(srtext);
	}
	else
	{
		elog(ERROR, "Cannot find SRID (%d) in spatial_ref_sys", srid);
	}

	spi_result = SPI_finish();
	if (spi_result != SPI_OK_FINISH)
	{
		elog(ERROR, "Could not disconnect from database using SPI");
	}

	return strs;
}


/**
 *  Given an SRID, return the corresponding proj strings
 *  (auth_name:auth_srid/srtext/proj4text)
 *
 *  If the integer is one of the "well known" projections we support
 *  (WGS84 UTM N/S, Polar Stereographic N/S - see SRID_* macros),
 *  return the proj4text for those.
 */
static PjStrs
GetProjStrings(int32_t srid)
{
	PjStrs strs;
	memset(&strs, 0, sizeof(strs));

	/* SRIDs in SPATIAL_REF_SYS */
	if ( srid < SRID_RESERVE_OFFSET )
	{
		return GetProjStringsSPI(srid);
	}
	/* Automagic SRIDs */
	else
	{
		strs.proj4text = palloc(maxprojlen);
		int id = srid;
		/* UTM North */
		if ( id >= SRID_NORTH_UTM_START && id <= SRID_NORTH_UTM_END )
		{
			snprintf(strs.proj4text, maxprojlen, "+proj=utm +zone=%d +ellps=WGS84 +datum=WGS84 +units=m +no_defs", id - SRID_NORTH_UTM_START + 1);
		}
		/* UTM South */
		else if ( id >= SRID_SOUTH_UTM_START && id <= SRID_SOUTH_UTM_END )
		{
			snprintf(strs.proj4text, maxprojlen, "+proj=utm +zone=%d +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs", id - SRID_SOUTH_UTM_START + 1);
		}
		/* Lambert zones (about 30x30, larger in higher latitudes) */
		/* There are three latitude zones, divided at -90,-60,-30,0,30,60,90. */
		/* In yzones 2,3 (equator) zones, the longitudinal zones are divided every 30 degrees (12 of them) */
		/* In yzones 1,4 (temperate) zones, the longitudinal zones are every 45 degrees (8 of them) */
		/* In yzones 0,5 (polar) zones, the longitudinal zones are ever 90 degrees (4 of them) */
		else if ( id >= SRID_LAEA_START && id <= SRID_LAEA_END )
		{
			int zone = id - SRID_LAEA_START;
			int xzone = zone % 20;
			int yzone = zone / 20;
			double lat_0 = 30.0 * (yzone - 3) + 15.0;
			double lon_0 = 0.0;

			/* The number of xzones is variable depending on yzone */
			if  ( yzone == 2 || yzone == 3 )
				lon_0 = 30.0 * (xzone - 6) + 15.0;
			else if ( yzone == 1 || yzone == 4 )
				lon_0 = 45.0 * (xzone - 4) + 22.5;
			else if ( yzone == 0 || yzone == 5 )
				lon_0 = 90.0 * (xzone - 2) + 45.0;
			else
				lwerror("Unknown yzone encountered!");

			snprintf(strs.proj4text, maxprojlen, "+proj=laea +ellps=WGS84 +datum=WGS84 +lat_0=%g +lon_0=%g +units=m +no_defs", lat_0, lon_0);
		}
		/* Lambert Azimuthal Equal Area South Pole */
		else if ( id == SRID_SOUTH_LAMBERT )
		{
			strncpy(strs.proj4text, "+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxprojlen );
		}
		/* Polar Sterographic South */
		else if ( id == SRID_SOUTH_STEREO )
		{
			strncpy(strs.proj4text, "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxprojlen);
		}
		/* Lambert Azimuthal Equal Area North Pole */
		else if ( id == SRID_NORTH_LAMBERT )
		{
			strncpy(strs.proj4text, "+proj=laea +lat_0=90 +lon_0=-40 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxprojlen );
		}
		/* Polar Stereographic North */
		else if ( id == SRID_NORTH_STEREO )
		{
			strncpy(strs.proj4text, "+proj=stere +lat_0=90 +lat_ts=71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxprojlen );
		}
		/* World Mercator */
		else if ( id == SRID_WORLD_MERCATOR )
		{
			strncpy(strs.proj4text, "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs", maxprojlen );
		}
		else
		{
			elog(ERROR, "Invalid reserved SRID (%d)", srid);
			return strs;
		}

		POSTGIS_DEBUGF(3, "returning on SRID=%d: %s", srid, strs.proj4text);
		return strs;
	}
}

static int
pjstrs_has_entry(const PjStrs *strs)
{
	if ((strs->proj4text && strlen(strs->proj4text)) ||
		(strs->authtext && strlen(strs->authtext)) ||
		(strs->srtext && strlen(strs->srtext)))
		return 1;
	else
		return 0;
}

static void
pjstrs_pfree(PjStrs *strs)
{
	if (strs->proj4text)
		pfree(strs->proj4text);
	if (strs->authtext)
		pfree(strs->authtext);
	if (strs->srtext)
		pfree(strs->srtext);
}

static char*
pgstrs_get_entry(const PjStrs *strs, int n)
{
	switch (n)
	{
		case 0:
			return strs->authtext;
		case 1:
			return strs->srtext;
		case 2:
			return strs->proj4text;
		default:
			return NULL;
	}
}

/**
 * Add an entry to the local PROJ SRS cache. If we need to wrap around then
 * we must make sure the entry we choose to delete does not contain other_srid
 * which is the definition for the other half of the transformation.
 */
static LWPROJ *
AddToPROJSRSCache(PROJSRSCache *PROJCache, int32_t srid_from, int32_t srid_to)
{
	MemoryContext oldContext;

	PjStrs from_strs, to_strs;
	char *pj_from_str, *pj_to_str;

	/*
	** Turn the SRID number into a proj4 string, by reading from spatial_ref_sys
	** or instantiating a magical value from a negative srid.
	*/
	from_strs = GetProjStrings(srid_from);
	if (!pjstrs_has_entry(&from_strs))
		elog(ERROR, "got NULL for SRID (%d)", srid_from);
	to_strs = GetProjStrings(srid_to);
	if (!pjstrs_has_entry(&to_strs))
		elog(ERROR, "got NULL for SRID (%d)", srid_to);

	oldContext = MemoryContextSwitchTo(PROJCache->PROJSRSCacheContext);


	LWPROJ *projection = NULL;
	/* Try combinations of AUTH_NAME:AUTH_SRID/SRTEXT/PROJ4TEXT until we find */
	/* one that gives us a usable transform. Note that we prefer */
	/* AUTH_NAME:AUTH_SRID over SRTEXT and SRTEXT over PROJ4TEXT */
	/* (3 entries * 3 entries = 9 combos) */
	uint32_t i;
	for (i = 0; i < 9; i++)
	{
		pj_from_str = pgstrs_get_entry(&from_strs, i / 3);
		pj_to_str   = pgstrs_get_entry(&to_strs,   i % 3);
		if (!(pj_from_str && pj_to_str))
			continue;

		projection = lwproj_from_str(pj_from_str, pj_to_str);
		if (projection)
			break;
	}
	if (!projection)
	{
		elog(ERROR, "could not form projection (LWPROJ) from 'srid=%d' to 'srid=%d'", srid_from, srid_to);
		return NULL;
	}

	/* If the cache is already full then find the least used element and delete it */
	uint32_t cache_position = PROJCache->PROJSRSCacheCount;
	uint32_t hits = 1;
	if (cache_position == PROJ_CACHE_ITEMS)
	{
		cache_position = 0;
		hits = PROJCache->PROJSRSCache[0].hits;
		for (uint32_t i = 1; i < PROJ_CACHE_ITEMS; i++)
		{
			if (PROJCache->PROJSRSCache[i].hits < hits)
			{
				cache_position = i;
				hits = PROJCache->PROJSRSCache[i].hits;
			}
		}
		DeleteFromPROJSRSCache(PROJCache, cache_position);
		/* To avoid the element we are introduced now being evicted next (as
		 * it would have 1 hit, being most likely the lower one) we reuse the
		 * hits from the evicted position and add some extra buffer
		 */
		hits += 5;
	}
	else
	{
		PROJCache->PROJSRSCacheCount++;
	}

	POSTGIS_DEBUGF(3,
		       "adding transform %d => %d aka \"%s\" => \"%s\" to query cache at index %d",
		       srid_from,
		       srid_to,
		       pj_from_str,
		       pj_to_str,
		       cache_position);

	/* Free the projection strings */
	pjstrs_pfree(&from_strs);
	pjstrs_pfree(&to_strs);

	/* Store everything in new cache entry */
	PROJCache->PROJSRSCache[cache_position].srid_from = srid_from;
	PROJCache->PROJSRSCache[cache_position].srid_to = srid_to;
	PROJCache->PROJSRSCache[cache_position].projection = projection;
	PROJCache->PROJSRSCache[cache_position].hits = hits;

	MemoryContextSwitchTo(oldContext);
	return projection;
}

static void
DeleteFromPROJSRSCache(PROJSRSCache *PROJCache, uint32_t position)
{
	POSTGIS_DEBUGF(3,
		       "removing query cache entry with SRIDs %d => %d at index %u",
		       PROJCache->PROJSRSCache[position].srid_from,
		       PROJCache->PROJSRSCache[position].srid_to,
		       position);

	/* Call PROJSRSDestroyPJ to free the PROJ objects memory now instead of
	 * waiting for the parent memory context to exit */
	PROJSRSDestroyPJ(PROJCache->PROJSRSCache[position].projection);

	PROJCache->PROJSRSCache[position].projection = NULL;
	PROJCache->PROJSRSCache[position].srid_from = SRID_UNKNOWN;
	PROJCache->PROJSRSCache[position].srid_to = SRID_UNKNOWN;
}


int
lwproj_lookup(int32_t srid_from, int32_t srid_to, LWPROJ **pj)
{
	/* get or initialize the cache for this round */
	PROJSRSCache* proj_cache = GetPROJSRSCache();
	if (!proj_cache)
		return LW_FAILURE;

	postgis_initialize_cache();
	/* Add the output srid to the cache if it's not already there */
	*pj = GetProjectionFromPROJCache(proj_cache, srid_from, srid_to);
	if (*pj == NULL)
	{
		*pj = AddToPROJSRSCache(proj_cache, srid_from, srid_to);
	}

	return *pj != NULL;
}

int
lwproj_is_latlong(const LWPROJ *pj)
{
	return pj->source_is_latlong;
}

static int
srid_is_latlong(int32_t srid)
{
	LWPROJ *pj;
	if ( lwproj_lookup(srid, srid, &pj) == LW_FAILURE)
		return LW_FALSE;
	return lwproj_is_latlong(pj);
}

void
srid_check_latlong(int32_t srid)
{
	if (srid == SRID_DEFAULT || srid == SRID_UNKNOWN)
		return;

	if (srid_is_latlong(srid))
		return;

	ereport(ERROR, (
	            errcode(ERRCODE_INVALID_PARAMETER_VALUE),
	            errmsg("Only lon/lat coordinate systems are supported in geography.")));
}

srs_precision
srid_axis_precision(int32_t srid, int precision)
{
	srs_precision sp;
	sp.precision_xy = precision;
	sp.precision_z = precision;
	sp.precision_m = precision;

	if ( srid == SRID_UNKNOWN )
		return sp;

	if ( srid_is_latlong(srid) )
	{
		sp.precision_xy += 5;
		return sp;
	}

	return sp;
}

int
spheroid_init_from_srid(int32_t srid, SPHEROID *s)
{
	LWPROJ *pj;

	if ( lwproj_lookup(srid, srid, &pj) == LW_FAILURE)
		return LW_FAILURE;

	if (!pj->source_is_latlong)
		return LW_FAILURE;
	spheroid_init(s, pj->source_semi_major_metre, pj->source_semi_minor_metre);

	return LW_SUCCESS;
}