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;
}

static inline int
get_u16 (const U8 *ptr)
{
  return ptr[0] | (ptr[1] << 8);
}

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 = get_u16 (ptr + 0) * ( 90. / 32767.);
                    double plon = get_u16 (ptr + 2) * (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);
        sv_grow (self, sizeof (struct cdb));
        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.5
Committed:	Thu Mar 17 00:32:54 2022 UTC (3 years, 11 months ago) by root
Branch:	MAIN
CVS Tags:	rel-0_9
Changes since 1.4:	+11 -4 lines
Log Message:	0.9
#	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.5	static inline int
39			get_u16 (const U8 *ptr)
40			{
41			return ptr[0] \| (ptr[1] << 8);
42			}
43
44	root	1.1	MODULE = Geo::LatLon2Place PACKAGE = Geo::LatLon2Place
45
46			PROTOTYPES: ENABLE
47
48	root	1.4	SV *
49			lookup_ext_ (SV *cdb, int km, int boxes, NV lat, NV lon, int r0, int r1, int flags)
50	root	1.1	CODE:
51	root	1.4	{
52			struct cdb db = (struct cdb )SvPVX (cdb);
53
54			if (!r1)
55			r1 = km;
56
57			r0 = r0 / km;
58			r1 = (r1 + km - 1) / km;
59			double coslat = cos (TORAD (lat));
60			int cy = (lat + 90.) * boxes * (1. / 180.);
61			int x, y;
62
63			if (r1 < r0 \|\| r0 < 0 \|\| r1 < 0 \|\| r0 >= boxes / 2 \|\| r1 >= boxes / 2)
64			XSRETURN_EMPTY;
65
66			if (lat < -90. \|\| lat > 90. \|\| lon < -180 \|\| lon > 180.)
67			XSRETURN_EMPTY;
68
69			double mind = 1e99;
70			const U8 *resptr;
71			int reslen = 0;
72
73			for (y = intmax (0, cy - r1); y <= intmin (boxes - 1, cy + r1); ++y)
74			{
75			double glat = y * (180. / boxes) - 90.;
76			double coslat = cos (TORAD (glat));
77			int blat = boxes * coslat; /* can be zero */
78			int cx = (lon + 180.) * blat * (1. / 360.);
79
80			for (x = cx - r1; x <= cx + r1; ++x)
81			{
82			int rx = x;
83			rx += rx < 0 ? blat : 0;
84			rx -= rx >= blat ? blat : 0;
85
86			unsigned char key[4] = {
87			rx, rx >> 8,
88			y, y >> 8,
89			};
90
91	root	1.5	//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));
92	root	1.4
93			if (cdb_find (db, key, sizeof (key)) <= 0)
94			continue;
95
96			int len = cdb_datalen (db);
97			const U8 *ptr = cdb_get (db, len, cdb_datapos (db));
98
99			while (len > 0)
100			{
101			int datalen = ptr[5];
102
103	root	1.5	double plat = get_u16 (ptr + 0) * ( 90. / 32767.);
104			double plon = get_u16 (ptr + 2) * (180. / 32767.);
105	root	1.4	int w = ptr[4];
106
107			double dx = TORAD (lon - plon) * coslat;
108			double dy = TORAD (lat - plat);
109			double d2 = (dx * dx + dy * dy) * w;
110	root	1.5	//printf ("%g,%g %g %.*s\n", plat,plon,d2, datalen,ptr+6);
111	root	1.4
112			if (d2 < mind)
113			{
114			mind = d2;
115			resptr = ptr;
116			reslen = datalen;
117			}
118
119			len -= datalen + 6;
120			ptr += datalen + 6;
121			}
122			}
123			}
124
125			if (!reslen)
126			XSRETURN_EMPTY;
127	root	1.1
128	root	1.4	RETVAL = newSVpvn (resptr + 6, reslen);
129	root	1.1	}
130	root	1.4	OUTPUT: RETVAL
131	root	1.1
132	root	1.4	#############################################################################
133	root	1.1
134			int
135			cdb_init (SV *self, int fd)
136			CODE:
137			sv_upgrade (self, SVt_PV);
138	root	1.5	sv_grow (self, sizeof (struct cdb));
139	root	1.1	SvCUR_set (self, sizeof (struct cdb));
140			SvPOK_only (self);
141			RETVAL = cdb_init ((struct cdb *)SvPVX (self), fd);
142			OUTPUT: RETVAL
143
144			void
145			cdb_free (SV *self)
146			CODE:
147			cdb_free ((struct cdb *)SvPVX (self));
148
149			SV *
150			cdb_get (SV self, SV key)
151			CODE:
152			{
153			unsigned int p;
154			STRLEN l;
155			const char *s = SvPVbyte (key, l);
156			struct cdb db = (struct cdb )SvPVX (self);
157
158			if (cdb_find (db, s, l) <= 0)
159			XSRETURN_UNDEF;
160
161			p = cdb_datapos (db);
162			l = cdb_datalen (db);
163			RETVAL = newSVpvn (cdb_get (db, l, p), l);
164			}
165			OUTPUT: RETVAL
166
167	root	1.4	#############################################################################
168
169			int
170			cdb_make_start (int fd)
171			CODE:
172			RETVAL = cdb_make_start (&make, fd);
173			OUTPUT: RETVAL
174
175			int
176			cdb_make_add (SV k, SV v)
177			CODE:
178			{
179			STRLEN klen; const char *kp = SvPVbyte (k, klen);
180			STRLEN vlen; const char *vp = SvPVbyte (v, vlen);
181			RETVAL = cdb_make_add (&make, kp, klen, vp, vlen);
182			}
183			OUTPUT: RETVAL
184
185			int
186			cdb_make_finish ()
187			CODE:
188			RETVAL = cdb_make_finish (&make);
189			OUTPUT: RETVAL
190