/* Copyright (c) 2001, 2002, 2003 by Intevation GmbH
 * Authors:
 * Bernhard Herzog <bh@intevation.de>
 *
 * This program is free software under the GPL (>=v2)
 * Read the file COPYING coming with Thuban for details.
 */

/* module to read a shape from a shapefile with shapelib, project the
 * data with proj4 and draw it on a wxWindows DC
 */

#include <stdlib.h>
#include <math.h>

#include <Python.h>

#include <shapefil.h>
#include "pyshapelib_api.h"

/* PVALUE is a typef in projects.h and is used in winreg.h on Windows as
 * well. We don't need it, so we temporarily rename the one from
 * projects.h */
#define PVALUE PROJ_PVALUE

#include <projects.h>

/* undefine it again */
#undef PVALUE

#include <wx/wx.h>


static PyShapeLibAPI * pyshapelib_api;

/* helper function to extract the pointer value from a swig ptr string
 */
static void *
decode_pointer(char *string)
{
  unsigned long p = 0;

  /* Pointer values must start with leading underscore */
  if (*string == '_')
  {
      string++;
      /* Extract hex value from pointer */
      while (*string)
      {
	  if ((*string >= '0') && (*string <= '9'))
	    p = (p << 4) + (*string - '0');
	  else if ((*string >= 'a') && (*string <= 'f'))
	    p = (p << 4) + ((*string - 'a') + 10);
	  else
	    break;
	  string++;
      }
  }
  return (void*)p;
}

/* Return the object wrapped by the SWIG shadow object shadow */
static void*
get_pointer(PyObject* shadow)
{
    PyObject * string;
    void * result = NULL;

    string = PyObject_GetAttrString(shadow, "this");
    if (string && PyString_Check(string))
    {
	result = decode_pointer(PyString_AsString(string));
    }
    Py_XDECREF(string);

    return result;
}

/* Extract a pointer from a python object that has a cobject method.
 *
 * Call the object's cobject method and if it returns a Py_CObject
 * extract the pointer from it and write the it into *output and return
 * true. If unsuccessful set a Python exception, do not modify *output
 * and return false.
 *
 * If object is None, set *output will to NULL.
 */
template<class T>
int
extract_pointer(PyObject * object, T**output)
{
    void * value = NULL;
    PyObject *cobject = NULL;
    
    if (Py_None == object)
    {
	value = NULL;
    }
    else
    {
	cobject = PyObject_CallMethod(object, "cobject", NULL);
	if (!cobject)
	    return 0;

	/* now cobject should be a CObject containing the PJ* pointer */
	if (PyCObject_Check(cobject))
	{
	    value = PyCObject_AsVoidPtr(cobject);
	}
	else
	{
	    PyErr_SetString(PyExc_TypeError,
			    "The projection must be either None, "
			    "or an object whose cobject method "
			    "returns a CObjecta cobject");
	    Py_XDECREF(cobject);
	    return 0;
	}
    }

    *output = (T*)value;
    return 1;
}


/* Project the point (x, y) with the projections forward and inverse and
 * scale it with scalex, scaley and offset it by offx, offy and store in
 * *xdest, and *ydest.
 */
static void
project_point(double * xdest, double *ydest,
	      PJ* forward, PJ* inverse,
	      double scalex, double scaley, double offx, double offy,
	      double x, double y)
{
    projUV proj_point;

    proj_point.u = x; 
    proj_point.v = y;

    if (inverse)
        proj_point = pj_inv(proj_point, inverse);

    if (forward)
    {
	if (!inverse)
	{
	    /* If a map projection is used but no layer projection we
	     * have to assume that the shapefille data is in geographic
	     * coordinates. Thuban assumes that they're in degrees but
	     * proj uses radians, so we have to convert explicitly. With
	     * a layer projection, i.e. inverse, present we already have
	     * the coordinates in radians.
	     */
	    proj_point.u *= DEG_TO_RAD;
	    proj_point.v *= DEG_TO_RAD;
	}
        proj_point = pj_fwd(proj_point, forward);
    }
    else if (inverse)
    {
	/* if there's an inverse projection but no forward projection
	 * the coordinates are now in radians, so we have to convert
	 * them to degrees.
	 */
	proj_point.u *= RAD_TO_DEG;
	proj_point.v *= RAD_TO_DEG;
    }
	

    *xdest = (scalex * proj_point.u + offx);
    *ydest = (scaley * proj_point.v + offy);
}


/* Return the projected points as a wxPoint array. The array returned
 * has a length of num_vertices + num_parts - 1
 *
 * The return value is suitable as a parameter to DrawPolygon with a
 * non-null brush so that the shape will be filled correctly, even
 * shapes with holes or otherwise consisting of multiple parts. This is
 * achieved by connecting the start/end points of the parts with lines
 * in such a way that these lines are given twice, but in opposite
 * directions. This trick works on both X11 and Windows 2000 (perhaps
 * on other windows versions too). This trick is the reason why the
 * returned array has num_vertices + num_parts - 1 items.
 */
static wxPoint *
project_points(int num_vertices, int num_parts,
	       double * xs, double * ys, int * part_start,
	       PJ* forward, PJ* inverse,
	       double scalex, double scaley, double offx, double offy)
{
    int i;
    int num_points = num_vertices + num_parts - 1;
    wxPoint* points = (wxPoint*)malloc(num_points * sizeof(wxPoint));
    if (!points)
    {
	PyErr_NoMemory();
	return NULL;
    }

    for (i = 0; i < num_vertices; i++)
    {
        double x, y;

        project_point(&x, &y, forward, inverse,
                      scalex, scaley, offx, offy,
                      xs[i], ys[i]);
        points[i].x = (int)x;
        points[i].y = (int)y;
    }

    /* link all the parts so that we can draw a multipart polygon with
     * one DrawPolgon call. If the points array as a whole is used for
     * DrawPolgon, the connections between the start points of the
     * subparts are already in the polygon, but only in one direction.
     * Here we add the same connections, but in the opposite order and
     * direction.
     */
    for (i = num_parts - 1; i > 0; i--)
    {
	points[num_vertices + num_parts - 1 - i] = points[part_start[i]];
    }

    return points;
}

/*
 * Structure to hold common data for calls to draw_polygon_shape().
 */
struct s_draw_info
{
    PyObject * py_shapefile;
    PyObject * py_forward;
    PyObject * py_inverse;
    PyObject * py_dc;
    double scalex, scaley, offx, offy;
    SHPHandle handle;
    PJ * forward;
    PJ * inverse;
    wxDC * dc;
};


/* Free a draw info struct */
static void
free_draw_info(void * p_draw_info)
{
    s_draw_info * draw_info = (s_draw_info*)p_draw_info;
    Py_XDECREF(draw_info->py_shapefile);
    Py_XDECREF(draw_info->py_dc);
    Py_XDECREF(draw_info->py_forward);
    Py_XDECREF(draw_info->py_inverse);
    delete draw_info;
}


/*
 * Initialize the structure for drawing polygons. Call once for each
 * layer and pass the return value to draw_polygon_shape.
 */
static PyObject *
draw_polygon_init(PyObject *, PyObject * args)
{
    PyObject * py_shapefile;
    PyObject * py_forward;
    PyObject * py_inverse;
    PyObject * py_dc;
    double scalex, scaley, offx, offy;

    if (!PyArg_ParseTuple(args, "OOOOdddd",
			  &py_shapefile,
                          &py_dc, &py_forward, &py_inverse,
			  &scalex, &scaley, &offx, &offy))
	return NULL;

    s_draw_info * draw_info = new s_draw_info;

    draw_info->py_shapefile = 0;
    draw_info->py_dc = 0;
    draw_info->py_forward = 0;
    draw_info->py_inverse = 0;

    draw_info->py_shapefile = py_shapefile;
    Py_INCREF(draw_info->py_shapefile);
    if (!extract_pointer(draw_info->py_shapefile, &draw_info->handle))
	goto fail;

    draw_info->py_dc = py_dc;
    Py_INCREF(draw_info->py_dc);
    draw_info->dc = (wxDC*)get_pointer(draw_info->py_dc);
    if (!draw_info->dc)
    {
	PyErr_SetString(PyExc_TypeError,
			"third argument must be a wxDC instance");
	goto fail;
    }

    draw_info->py_forward = py_forward;
    Py_INCREF(draw_info->py_forward);
    if (!extract_pointer(draw_info->py_forward, &draw_info->forward))
	goto fail;

    draw_info->py_inverse = py_inverse;
    Py_INCREF(draw_info->py_inverse);
    if (!extract_pointer(draw_info->py_inverse, &draw_info->inverse))
	goto fail;

    draw_info->scalex = scalex;
    draw_info->scaley = scaley;
    draw_info->offx = offx;
    draw_info->offy = offy;

    return PyCObject_FromVoidPtr(draw_info, free_draw_info);

 fail:
    free_draw_info(draw_info);
    return NULL;
}


/*
 * Draw a polygon. draw_polygon_init() MUST be called prior to this function.
 */
static PyObject *
draw_polygon_shape(PyObject *, PyObject * args)
{
    int shape_index = 0;
    PyObject * brush_shadow;
    PyObject * pen_shadow;
    SHPObject * shape;
    wxPoint* points;
    wxPen * pen;
    wxBrush * brush;
    int num_points;
    wxDC * dc;
    PyObject * py_draw_info;
    s_draw_info * draw_info;

    if (!PyArg_ParseTuple(args, "O!iOO",
			  &PyCObject_Type, &py_draw_info,
			  &shape_index, &pen_shadow, &brush_shadow))
	return NULL;

    draw_info = (s_draw_info*)PyCObject_AsVoidPtr(py_draw_info);
    
    if (pen_shadow == Py_None)
    {
	pen = NULL;
    }
    else
    {
	pen = (wxPen*)get_pointer(pen_shadow);
	if (!pen)
	{
	    PyErr_SetString(PyExc_TypeError,
		       "fourth argument must be a wxPen instance or None");
	    return NULL;
	}
    }

    if (brush_shadow == Py_None)
    {
	brush = NULL;
    }
    else
    {
	brush = (wxBrush*)get_pointer(brush_shadow);
	if (!brush)
	{
	    PyErr_SetString(PyExc_TypeError,
			"fifth argument must be a wxBrush instance or None");
	    return NULL;
	}
    }

    shape = pyshapelib_api->SHPReadObject(draw_info->handle, shape_index);
    if (!shape)
	return PyErr_Format(PyExc_ValueError,
			    "Can't get shape %d from shapefile", shape_index);
    
    num_points = shape->nVertices + shape->nParts - 1;
    points = project_points(shape->nVertices, shape->nParts,
			    shape->padfX, shape->padfY, shape->panPartStart,
			    draw_info->forward, 
                            draw_info->inverse, 
                            draw_info->scalex, 
                            draw_info->scaley, 
                            draw_info->offx, 
                            draw_info->offy);
    dc = draw_info->dc;

    // If the shape is a polygon and a non-transparent brush was given,
    // fill the polygon
    if (shape->nSHPType == SHPT_POLYGON
	&& brush && brush != wxTRANSPARENT_BRUSH)
    {
	dc->SetPen(*wxTRANSPARENT_PEN);
	dc->SetBrush(*brush);
	dc->DrawPolygon(num_points, points, 0, 0);
    }

    // If a non-transparent pen was given, stroke the outline
    if (pen && pen != wxTRANSPARENT_PEN)
    {
	draw_info->dc->SetPen(*pen);
	draw_info->dc->SetBrush(*wxTRANSPARENT_BRUSH);
	for (int i = 0; i < shape->nParts; i++)
	{
	    int length;

	    if (i < shape->nParts - 1)
		length = shape->panPartStart[i + 1] - shape->panPartStart[i];
	    else
		length = shape->nVertices - shape->panPartStart[i];

	    // If the shape is a polygon, use DrawPolygon else use DrawLines
	    if (shape->nSHPType == SHPT_POLYGON)
	    {
		dc->DrawPolygon(length, points + shape->panPartStart[i], 0, 0);
	    }
	    else
	    {
		dc->DrawLines(length, points + shape->panPartStart[i], 0, 0);
	    }
	}
    }
    
    free(points);

    pyshapelib_api->SHPDestroyObject(shape);

    Py_INCREF(Py_None);
    return Py_None;
}



/* determine whether the line fom (SX, SY) to (EX, EY) is `hit' by a
 * click at (PX, PY), or whether a polygon containing this line is hit
 * in the interior at (PX, PY).
 *
 * Return -1 if the line it self his hit. Otherwise, return +1 if a
 * horizontal line from (PX, PY) to (-Infinity, PY) intersects the line
 * and 0 if it doesn't.
 *
 * The nonnegative return values can be used to determine whether (PX, PY)
 * is an interior point of a polygon according to the even-odd rule.
 */

#define PREC_BITS 4

static int
test_line(int sx, int sy, int ex, int ey, int px, int py)
{
    long vx, vy, dx, dy, len, dist, not_horizontal;

    if (ey < sy)
    {
	dist = ex; ex = sx; sx = dist;
	dist = ey; ey = sy; sy = dist;
    }
    not_horizontal = ey > sy + (2 << PREC_BITS);
    if (not_horizontal && (py >= ey || py < sy))
	return 0;

    vx = ex - sx; vy = ey - sy;

    len = (long)sqrt(vx * vx + vy * vy);
    if (!len)
	/* degenerate case of coincident end points. Assumes that some
	 * other part of the code has already determined whether the end
	 * point is hit.
	 */
	return 0;

    dx = px - sx; dy = py - sy;
    dist = vx * dy - vy * dx;
    if ((not_horizontal || (px >= sx && px <= ex) || (px >= ex && px <= sx))
	&& abs(dist) <= (len << (PREC_BITS + 1)))
	return -1;

    /* horizontal lines (vy == 0) always return 0 here. */
    return vy && py < ey && py >= sy && dx * abs(vy) > vx * abs(dy);
}



static PyObject *
point_in_polygon_shape(PyObject *, PyObject * args)
{
    int shape_index = 0;
    PyObject * shphandle_cobject;
    PyObject * forward_cobject;
    PyObject * inverse_cobject;
    int filled, stroked;
    double scalex, scaley, offx, offy;
    int px, py;
    SHPHandle handle;
    PJ * forward = NULL;
    PJ * inverse = NULL;
    SHPObject * shape;
    wxPoint* points;
    int num_points;
    int cross_count, linehit;
    int result;

    if (!PyArg_ParseTuple(args, "O!iiiOOddddii",
			  &PyCObject_Type, &shphandle_cobject,
			  &shape_index, &filled, &stroked,
			  &forward_cobject, &inverse_cobject,
			  &scalex, &scaley, &offx, &offy, &px, &py))
	return NULL;

    handle = (SHPHandle)PyCObject_AsVoidPtr(shphandle_cobject);

    if (!extract_pointer(forward_cobject, &forward))
	return NULL;

    if (!extract_pointer(inverse_cobject, &inverse))
	return NULL;

    shape = pyshapelib_api->SHPReadObject(handle, shape_index);
    if (!shape)
	return PyErr_Format(PyExc_ValueError,
			    "Can't get shape %d from shapefile", shape_index);

    num_points = shape->nVertices + shape->nParts - 1;
    points = project_points(shape->nVertices, shape->nParts,
			    shape->padfX, shape->padfY, shape->panPartStart,
			    forward, inverse, scalex, scaley, offx, offy);

    long scaled_px = (px << PREC_BITS) + 1, scaled_py = (py << PREC_BITS) + 1;

    cross_count = 0; linehit = 0;
    for (int part = 0; part < shape->nParts; part++)
    {
	int start, end;

	if (part < shape->nParts - 1)
	{
	    start = shape->panPartStart[part];
	    end = shape->panPartStart[part + 1] - 1;
	}
	else
	{
	    start = shape->panPartStart[part];
	    end = shape->nVertices - 1;
	}

	for (int vertex = start; vertex < end; vertex++)
	{
	    //printf("test_line: %d, %d: %d, %d -- %d, %d",
	    //	   px, py, points[vertex].x, points[vertex].y,
	    //	   points[vertex + 1].x, points[vertex + 1].y);
	    result = test_line(points[vertex].x << PREC_BITS,
				points[vertex].y << PREC_BITS,
				points[vertex + 1].x << PREC_BITS,
				points[vertex + 1].y << PREC_BITS,
				scaled_px, scaled_py);
	    //printf(":\t%d\n", result);
	    if (result < 0)
	    {
		linehit = 1;
		break;
	    }
	    cross_count += result;
	}
	if (linehit)
	    break;
    }

    free(points);

    pyshapelib_api->SHPDestroyObject(shape);

    if (filled)
    {
	if (stroked && linehit)
	    result = -1;
	else
	    result = cross_count % 2;
    }
    else if (stroked)
    {
	if (linehit)
	    result = -1;
	else
	    result = 0;
    }
    else
	result = 0;
		

    return PyInt_FromLong(result);
}


static PyObject*
shape_centroid(PyObject *self, PyObject *args)
{
    int shape_index = 0;
    PyObject * shphandle_cobject;
    PyObject * forward_cobject;
    PyObject * inverse_cobject;
    double scalex, scaley, offx, offy;
    SHPHandle handle;
    PJ * forward = NULL;
    PJ * inverse = NULL;
    SHPObject * shape;
    double centroidx = 0, centroidy = 0;
    double sum = 0, area;
    double lastx=0, lasty=0, x, y;

    if (!PyArg_ParseTuple(args, "O!iOOdddd",
			  &PyCObject_Type, &shphandle_cobject,
			  &shape_index, &forward_cobject, &inverse_cobject,
			  &scalex, &scaley, &offx, &offy))
	return NULL;

    handle = (SHPHandle)PyCObject_AsVoidPtr(shphandle_cobject);

    if (!extract_pointer(forward_cobject, &forward))
	return NULL;

    if (!extract_pointer(inverse_cobject, &inverse))
	return NULL;

    shape = pyshapelib_api->SHPReadObject(handle, shape_index);
    if (!shape)
	return PyErr_Format(PyExc_ValueError,
			    "Can't get shape %d from shapefile", shape_index);

    for (int part = 0; part < shape->nParts; part++)
    {
	int start, end;

	if (part < shape->nParts - 1)
	{
	    start = shape->panPartStart[part];
	    end = shape->panPartStart[part + 1];
	}
	else
	{
	    start = shape->panPartStart[part];
	    end = shape->nVertices;
	}

	project_point(&lastx, &lasty, forward, inverse,
		      scalex, scaley, offx, offy,
		      shape->padfX[start],
		      shape->padfY[start]);
	for (int vertex = start + 1; vertex < end; vertex++)
	{
	    project_point(&x, &y, forward, inverse,
			  scalex, scaley, offx, offy,
			  shape->padfX[vertex],
			  shape->padfY[vertex]);
	    area = x * lasty - lastx * y;
	    centroidx += area * (x + lastx);
	    centroidy += area * (y + lasty);
	    sum += area;
	    lastx = x;
	    lasty = y;
	}
    }

    pyshapelib_api->SHPDestroyObject(shape);

    // Handle pathologic case where all points are equal.
    if (sum == 0.0)
    	return Py_BuildValue("dd", lastx, lasty);
    else
    	return Py_BuildValue("dd", centroidx / (3 * sum), centroidy / (3 * sum));
}


/*
 * Returns a tuple containing the version of PROJ that was used when
 * compiling this file (major, minor, release). If the version was
 * not available (because of a very old version of PROJ) the tuple
 * will be empty.
 */

static PyObject*
get_proj_version(PyObject *, PyObject * args)
{
    PyObject * result;

#ifdef PJ_VERSION
    result = PyTuple_New(3);
    if (result == NULL)
        return NULL;

    int ver = PJ_VERSION;
    int major = ver / 100; ver -= major * 100;
    int minor = ver / 10;  ver -= minor * 10;
    int release = ver;

    PyTuple_SET_ITEM(result, 0, PyInt_FromLong(major));
    PyTuple_SET_ITEM(result, 1, PyInt_FromLong(minor));
    PyTuple_SET_ITEM(result, 2, PyInt_FromLong(release));
#else
    result = PyTuple_New(0);
    if (result == NULL)
        return NULL;
#endif

    Py_INCREF(result);
    return result;
}


static PyObject*
get_gtk_version(PyObject *, PyObject * args)
{
    PyObject * result;

#ifdef __WXGTK__
    extern const int gtk_major_version;
    extern const int gtk_minor_version;
    extern const int gtk_micro_version;

    result = PyTuple_New(3);
    if (result == NULL)
        return NULL;

    PyTuple_SET_ITEM(result, 0, PyInt_FromLong(gtk_major_version));
    PyTuple_SET_ITEM(result, 1, PyInt_FromLong(gtk_minor_version));
    PyTuple_SET_ITEM(result, 2, PyInt_FromLong(gtk_micro_version));
#else
    result = PyTuple_New(0);
    if (result == NULL)
        return NULL;
#endif

    Py_INCREF(result);
    return result;
}

static PyMethodDef wxproj_methods[] = {
    { "draw_polygon_shape", draw_polygon_shape, METH_VARARGS },
    { "draw_polygon_init", draw_polygon_init, METH_VARARGS },
    { "point_in_polygon_shape", point_in_polygon_shape, METH_VARARGS },
    { "shape_centroid", shape_centroid, METH_VARARGS },
    { "get_proj_version", get_proj_version, METH_VARARGS },
    { "get_gtk_version", get_gtk_version, METH_VARARGS },
    {NULL, NULL}
};


#ifdef __cplusplus
extern "C" 
#endif
void initwxproj(void)
{
    PyObject * shapelib = NULL;
    PyObject * c_api_func = NULL;
    PyObject * cobj = NULL;

    Py_InitModule("wxproj", wxproj_methods);

    PYSHAPELIB_IMPORT_API(pyshapelib_api);

    Py_XDECREF(cobj);
    Py_XDECREF(c_api_func);
    Py_XDECREF(shapelib);
}



