Geo-LatLon2Place/LatLon2Place.xs

#include "EXTERN.h"
#include "perl.h"
#include "XSUB.h"

#include <math.h>

#include "perlmulticore.h"

#if EMBED_CDB
 #include "cdb-embedded.c"
#else
 #include <cdb.h>
#endif

#define TORAD(r) ((r) * (M_PI / 180.))

static struct cdb_make make;

struct res
{
  double mind;
  unsigned int respos, reslen;
  double x, y;
};

static inline int
intmin (int a, int b)
{
  return a > b ? b : a;
}

static inline int
intmax (int a, int b)
{
  return a > b ? a : b;
}

MODULE = Geo::LatLon2Place              PACKAGE = Geo::LatLon2Place

PROTOTYPES: ENABLE

SV *
lookup_ext_ (SV *cdb, int km, int boxes, NV lat, NV lon, int r0, int r1, int flags)
        CODE:
{
        struct cdb *db = (struct cdb *)SvPVX (cdb);

        if (!r1)
          r1 = km;

        r0 =  r0           / km;
        r1 = (r1 + km - 1) / km;
        double coslat = cos (TORAD (lat));
        int cy = (lat + 90.) * boxes * (1. / 180.);
        int x, y;

        if (r1 < r0 || r0 < 0 || r1 < 0 || r0 >= boxes / 2 || r1 >= boxes / 2)
          XSRETURN_EMPTY;

        if (lat < -90. || lat > 90. || lon < -180 || lon > 180.)
          XSRETURN_EMPTY;

        double mind = 1e99;
        const U8 *resptr;
        int reslen = 0;

        for (y = intmax (0, cy - r1); y <= intmin (boxes - 1, cy + r1); ++y)
          {
            double glat = y * (180. / boxes) - 90.;
            double coslat = cos (TORAD (glat));
            int blat = boxes * coslat; /* can be zero */
            int cx = (lon + 180.) * blat * (1. / 360.);

            for (x = cx - r1; x <= cx + r1; ++x)
              {
                int rx = x;
                rx += rx <     0 ? blat : 0;
                rx -= rx >= blat ? blat : 0;

                unsigned char key[4] = {
                  rx, rx >> 8,
                   y,  y >> 8,
                };

                printf ("x,y %4d,%4d blat %d %d %g %02x%02x%02x%02x %d\n", rx, y, blat, (int)glat, TORAD(glat), key[0],key[1],key[2],key[3], sizeof(key));

                if (cdb_find (db, key, sizeof (key)) <= 0)
                  continue;

                int len = cdb_datalen (db);
                const U8 *ptr = cdb_get (db, len, cdb_datapos (db));

                while (len > 0)
                  {
                    int datalen = ptr[5];

                    double plat = ptr[0] | (ptr[1] << 8) * ( 90. / 32767.);
                    double plon = ptr[2] | (ptr[3] << 8) * (180. / 32767.);
                    int w  = ptr[4];

                    double dx = TORAD (lon - plon) * coslat;
                    double dy = TORAD (lat - plat);
                    double d2 = (dx * dx + dy * dy) * w;
                    printf ("%g,%g %g %.*s\n", plat,plon,d2, datalen,ptr+6);

                    if (d2 < mind)
                      {
                        mind   = d2;
                        resptr = ptr;
                        reslen = datalen;
                      }

                    len -= datalen + 6;
                    ptr += datalen + 6;
                  }
              }
          }

        if (!reslen)
          XSRETURN_EMPTY;

        RETVAL = newSVpvn (resptr + 6, reslen);
}
        OUTPUT: RETVAL

#############################################################################

int
cdb_init (SV *self, int fd)
        CODE:
        sv_upgrade (self, SVt_PV);
        SvCUR_set (self, sizeof (struct cdb));
        SvPOK_only (self);
        RETVAL = cdb_init ((struct cdb *)SvPVX (self), fd);
        OUTPUT: RETVAL

void
cdb_free (SV *self)
        CODE:
        cdb_free ((struct cdb *)SvPVX (self));

SV *
cdb_get (SV *self, SV *key)
        CODE:
{
        unsigned int p;
        STRLEN l;
        const char *s = SvPVbyte (key, l);
        struct cdb *db = (struct cdb *)SvPVX (self);

        if (cdb_find (db, s, l) <= 0)
          XSRETURN_UNDEF;

        p = cdb_datapos (db);
        l = cdb_datalen (db);
        RETVAL = newSVpvn (cdb_get (db, l, p), l);
}
        OUTPUT: RETVAL

#############################################################################

int
cdb_make_start (int fd)
        CODE:
        RETVAL = cdb_make_start (&make, fd);
        OUTPUT: RETVAL

int
cdb_make_add (SV *k, SV *v)
        CODE:
{
        STRLEN klen; const char *kp = SvPVbyte (k, klen);
        STRLEN vlen; const char *vp = SvPVbyte (v, vlen);
        RETVAL = cdb_make_add (&make, kp, klen, vp, vlen);
}
        OUTPUT: RETVAL

int
cdb_make_finish ()
        CODE:
        RETVAL = cdb_make_finish (&make);
        OUTPUT: RETVAL

Revision:	1.4
Committed:	Tue Mar 15 07:33:40 2022 UTC (2 years, 2 months ago) by root
Branch:	MAIN
Changes since 1.3:	+128 -17 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	root	1.1	#include "EXTERN.h"
2			#include "perl.h"
3			#include "XSUB.h"
4
5	root	1.4	#include <math.h>
6
7	root	1.3	#include "perlmulticore.h"
8	root	1.1
9	root	1.3	#if EMBED_CDB
10			#include "cdb-embedded.c"
11			#else
12			#include <cdb.h>
13			#endif
14	root	1.2
15	root	1.4	#define TORAD(r) ((r) * (M_PI / 180.))
16
17	root	1.1	static struct cdb_make make;
18
19	root	1.4	struct res
20			{
21			double mind;
22			unsigned int respos, reslen;
23			double x, y;
24			};
25
26			static inline int
27			intmin (int a, int b)
28			{
29			return a > b ? b : a;
30			}
31
32			static inline int
33			intmax (int a, int b)
34			{
35			return a > b ? a : b;
36			}
37
38	root	1.1	MODULE = Geo::LatLon2Place PACKAGE = Geo::LatLon2Place
39
40			PROTOTYPES: ENABLE
41
42	root	1.4	SV *
43			lookup_ext_ (SV *cdb, int km, int boxes, NV lat, NV lon, int r0, int r1, int flags)
44	root	1.1	CODE:
45	root	1.4	{
46			struct cdb db = (struct cdb )SvPVX (cdb);
47
48			if (!r1)
49			r1 = km;
50
51			r0 = r0 / km;
52			r1 = (r1 + km - 1) / km;
53			double coslat = cos (TORAD (lat));
54			int cy = (lat + 90.) * boxes * (1. / 180.);
55			int x, y;
56
57			if (r1 < r0 \|\| r0 < 0 \|\| r1 < 0 \|\| r0 >= boxes / 2 \|\| r1 >= boxes / 2)
58			XSRETURN_EMPTY;
59
60			if (lat < -90. \|\| lat > 90. \|\| lon < -180 \|\| lon > 180.)
61			XSRETURN_EMPTY;
62
63			double mind = 1e99;
64			const U8 *resptr;
65			int reslen = 0;
66
67			for (y = intmax (0, cy - r1); y <= intmin (boxes - 1, cy + r1); ++y)
68			{
69			double glat = y * (180. / boxes) - 90.;
70			double coslat = cos (TORAD (glat));
71			int blat = boxes * coslat; /* can be zero */
72			int cx = (lon + 180.) * blat * (1. / 360.);
73
74			for (x = cx - r1; x <= cx + r1; ++x)
75			{
76			int rx = x;
77			rx += rx < 0 ? blat : 0;
78			rx -= rx >= blat ? blat : 0;
79
80			unsigned char key[4] = {
81			rx, rx >> 8,
82			y, y >> 8,
83			};
84
85			printf ("x,y %4d,%4d blat %d %d %g %02x%02x%02x%02x %d\n", rx, y, blat, (int)glat, TORAD(glat), key[0],key[1],key[2],key[3], sizeof(key));
86
87			if (cdb_find (db, key, sizeof (key)) <= 0)
88			continue;
89
90			int len = cdb_datalen (db);
91			const U8 *ptr = cdb_get (db, len, cdb_datapos (db));
92
93			while (len > 0)
94			{
95			int datalen = ptr[5];
96
97			double plat = ptr[0] \| (ptr[1] << 8) * ( 90. / 32767.);
98			double plon = ptr[2] \| (ptr[3] << 8) * (180. / 32767.);
99			int w = ptr[4];
100
101			double dx = TORAD (lon - plon) * coslat;
102			double dy = TORAD (lat - plat);
103			double d2 = (dx * dx + dy * dy) * w;
104			printf ("%g,%g %g %.*s\n", plat,plon,d2, datalen,ptr+6);
105
106			if (d2 < mind)
107			{
108			mind = d2;
109			resptr = ptr;
110			reslen = datalen;
111			}
112
113			len -= datalen + 6;
114			ptr += datalen + 6;
115			}
116			}
117			}
118
119			if (!reslen)
120			XSRETURN_EMPTY;
121	root	1.1
122	root	1.4	RETVAL = newSVpvn (resptr + 6, reslen);
123	root	1.1	}
124	root	1.4	OUTPUT: RETVAL
125	root	1.1
126	root	1.4	#############################################################################
127	root	1.1
128			int
129			cdb_init (SV *self, int fd)
130			CODE:
131			sv_upgrade (self, SVt_PV);
132			SvCUR_set (self, sizeof (struct cdb));
133			SvPOK_only (self);
134			RETVAL = cdb_init ((struct cdb *)SvPVX (self), fd);
135			OUTPUT: RETVAL
136
137			void
138			cdb_free (SV *self)
139			CODE:
140			cdb_free ((struct cdb *)SvPVX (self));
141
142			SV *
143			cdb_get (SV self, SV key)
144			CODE:
145			{
146			unsigned int p;
147			STRLEN l;
148			const char *s = SvPVbyte (key, l);
149			struct cdb db = (struct cdb )SvPVX (self);
150
151			if (cdb_find (db, s, l) <= 0)
152			XSRETURN_UNDEF;
153
154			p = cdb_datapos (db);
155			l = cdb_datalen (db);
156			RETVAL = newSVpvn (cdb_get (db, l, p), l);
157			}
158			OUTPUT: RETVAL
159
160	root	1.4	#############################################################################
161
162			int
163			cdb_make_start (int fd)
164			CODE:
165			RETVAL = cdb_make_start (&make, fd);
166			OUTPUT: RETVAL
167
168			int
169			cdb_make_add (SV k, SV v)
170			CODE:
171			{
172			STRLEN klen; const char *kp = SvPVbyte (k, klen);
173			STRLEN vlen; const char *vp = SvPVbyte (v, vlen);
174			RETVAL = cdb_make_add (&make, kp, klen, vp, vlen);
175			}
176			OUTPUT: RETVAL
177
178			int
179			cdb_make_finish ()
180			CODE:
181			RETVAL = cdb_make_finish (&make);
182			OUTPUT: RETVAL
183