--- /dev/null
+set role gis;
+
+drop schema if exists hex cascade;
+create schema hex;
+
+set search_path to hex,public,gis;
+
+create or replace function hexgrid(
+ xdim integer,
+ ydim integer,
+ angle numeric,
+ coverage geometry,
+ srid integer,
+ origin geometry(POINT),
+ tilesize numeric
+) returns table (x integer, y integer, cantor numeric, hex geometry) as $$
+declare
+ rcover geometry;
+ stride double precision;
+ x integer;
+ y integer;
+ xmin double precision;
+ ymin double precision;
+ xoff double precision;
+begin
+ -- the origin will be the center of the 1,1 hex.
+ -- if null, then the origin will be set to fit in the
+ -- covering box.
+
+ -- xdim, and ydim can be null, in which case they are
+ -- calculated from the tilesize and the covering box
+
+ -- calculate ydim
+
+ -- precedence: ydim, xdim, tilesize
+
+ -- counter rotate geometry, fill, rotate grid
+ rcover = ST_Rotate(st_transform(coverage,srid), -angle*3.1415926/180.0);
+
+ xmin := ST_XMin(rcover);
+ ymin := ST_YMin(rcover);
+
+ stride := tilesize * 2.0 / sqrt(3.0);
+
+ --raise notice 'stride = %', stride;
+ xdim := ((ST_XMax(rcover) - ST_XMin(rcover))/tilesize*2/sqrt(3.0))::integer;
+ ydim := ((ST_YMax(rcover) - ST_YMin(rcover))/tilesize)::integer;
+
+ --raise notice 'coverage = % % % %', xmin, ymin, xdim, ydim;
+
+ -- TODO reorient these so right hand rule is inside
+ hex := ST_Scale(ST_GeomFromText('POLYGON((
+.577350269189625764509148780502 0.0,
+.288675134594812882254574390251 0.5,
+ -.288675134594812882254574390251 0.5,
+ -.577350269189625764509148780502 0.0,
+ -.288675134594812882254574390251 -0.5,
+ .288675134594812882254574390251 -0.5,
+ .577350269189625764509148780502 0.0))',ST_SRID(rcover)),tilesize,tilesize,1.0);
+
+ xoff := .577350269189625764509148780502 * 1.5;
+
+ for x in 0 .. xdim loop
+ for y in 0 .. ydim loop
+ -- TODO actually calculate cantor
+ return query select x, ydim-y, hex.cantor(ARRAY[x,ydim-y]),
+ ST_Translate(hex,
+ tilesize * (x * xoff) + xmin,
+ --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
+ tilesize * (y + 0.5 * (x%2) ) + ymin
+ )
+ --st_rotate(ST_Translate(hex,
+ --tilesize * (x * xoff) + xmin,
+ --tilesize * (y + 0.5 * (x%2) - (x+1)/2) + ymin
+ --), -3.1415926/6)
+ ;
+ end loop;
+ end loop;
+END;
+$$
+language 'plpgsql'
+;
+
+-- create some hexagon ops
+
+create function ijk(xy integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ ijk integer[];
+begin
+ if array_length(xy,1) = 3 then
+ return xy;
+ end if;
+
+ ijk[1] := xy[1];
+ ijk[2] := -xy[2];
+
+ if xy[1] < 0 then
+ ijk[2] := ijk[2] + (-xy[1]+1) / 2;
+ else
+ ijk[2] := ijk[2] - xy[1]/2;
+ end if;
+
+ ijk[3] := -ijk[1] - ijk[2];
+
+ return ijk;
+end
+$$ immutable strict;
+
+create function ijk(x integer, y integer) returns integer[] language 'sql' as $$
+ select hex.ijk(ARRAY[x,y]);
+$$ immutable strict;
+
+-- invert the cantor paring
+create function xy(hex numeric) returns integer[] language 'plpgsql' as $$
+declare
+ xy integer[];
+ w integer;
+ t integer;
+ py integer;
+begin
+ if hex >= 0 then
+ hex := hex-1;
+ w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
+ t := (w * w + w)/2;
+ py := hex - t;
+
+ xy = ARRAY[w-py,py];
+ -- handle the negative cases
+ else
+ hex := -hex;
+ hex := hex-1;
+ w := floor((sqrt(8.0 * hex + 1.0) - 1.0)/2.0)::integer;
+ t := (w * w + w)/2;
+ py := hex - t;
+
+ xy = ARRAY[w-py,py];
+
+ if xy[1] % 2 = 1 then
+ xy[1] = -(xy[1]+1)/2;
+ else
+ xy[1] = xy[1]/2;
+ end if;
+ if xy[2] % 2 = 1 then
+ xy[2] = -(xy[2]+1)/2;
+ else
+ xy[2] = xy[2]/2;
+ end if;
+ end if;
+
+ return xy;
+end;
+$$ immutable strict;
+
+create function ijk(hex numeric) returns integer[] language 'sql' as $$
+ select hex.ijk(hex.xy(hex));
+$$ immutable strict;
+
+create function xy(ijk integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ xy integer[];
+begin
+ -- make sure this is really an ijk
+ -- just no-op if it's already an
+ -- xy, so this is sort of a coercion function
+ if array_length(ijk,1) = 2 then
+ return ijk;
+ end if;
+ xy[1] := ijk[1];
+ xy[2] := -ijk[2];
+
+ if ijk[1] < 0 then
+ xy[2] := xy[2] + (-ijk[1] + 1) / 2;
+ else
+ xy[2] := xy[2] - ijk[1]/2;
+ end if;
+ return xy;
+end;
+$$;
+
+CREATE OR REPLACE FUNCTION vec_add(arr1 integer[], arr2 integer[])
+RETURNS integer[] AS
+$$
+SELECT array_agg(result)
+FROM (SELECT tuple.val1 + tuple.val2 AS result
+ FROM (SELECT UNNEST($1) AS val1
+ ,UNNEST($2) AS val2
+ ,generate_subscripts($1, 1) AS ix) tuple
+ ORDER BY ix) inn;
+$$ LANGUAGE SQL IMMUTABLE STRICT;
+
+create function adjhex(dir integer, hex integer[]) returns integer[] language 'plpgsql' as $$
+declare
+ adj integer[];
+begin
+ if array_length(hex,1) = 2 then
+ hex = hex.ijk(hex);
+ end if;
+
+ dir := dir % 6;
+ if dir = 2 then
+ hex := hex.vec_add(hex, ARRAY[-1, 1, 0]);
+ elsif dir = 1 then
+ hex := hex.vec_add(hex, ARRAY[ 0, 1,-1]);
+ elsif dir = 0 then
+ hex := hex.vec_add(hex, ARRAY[ 1, 0,-1]);
+ elsif dir = 5 then
+ hex := hex.vec_add(hex, ARRAY[ 1,-1, 0]);
+ elsif dir = 4 then
+ hex := hex.vec_add(hex, ARRAY[ 0,-1, 1]);
+ elsif dir = 3 then
+ hex := hex.vec_add(hex, ARRAY[-1, 0, 1]);
+ else
+ -- error
+ end if;
+ return hex;
+end;
+$$ immutable strict;
+
+create function adjhex(dir integer, hex numeric) returns integer[] language 'sql' as $$
+ select hex.adjhex(dir, hex.ijk(hex));
+$$ immutable strict;
+
+-- TODO make this work for any, or inline it
+create function natcantor(hex integer[]) returns numeric language 'sql' as $$
+ select cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+$$ immutable strict;
+
+create function cantor(hex integer[]) returns numeric language 'plpgsql' as $$
+begin
+ -- convert to xy
+ if array_length(hex,1) = 3 then
+ hex = hex.xy(hex);
+ end if;
+
+ if hex[1] < 0 or hex[2] < 0 then
+ hex[1] = hex[1] * 2;
+ hex[2] = hex[2] * 2;
+ if hex[1] < 0 then
+ hex[1] = -hex[1] - 1;
+ end if;
+ if hex[2] < 0 then
+ hex[2] = -hex[2] - 1;
+ end if;
+ return -cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+ end if;
+
+ return cast((hex[1] + hex[2]) * (hex[1] + hex[2] + 1) / 2 + hex[2] + 1 as numeric);
+end;
+$$ immutable strict;
+
+create function adjacent(hex numeric) returns setof numeric language 'sql' as $$
+select hex.cantor(hex.adjhex(dir, hex)) from generate_series(0,5) as dir;
+$$ immutable strict;
+
+create function adjacent(dir integer, hex numeric) returns numeric language 'sql' as $$
+ select hex.cantor(hex.adjhex(dir, hex.ijk(hex)));
+$$ immutable strict;
+
+create function within(hex numeric, range integer = 1, include boolean = false) returns setof numeric language 'plpgsql' as $$
+declare
+ c integer[]; -- cantor coordinates
+ base integer[];
+begin
+ if range <= 0 then
+ return;
+ end if;
+ if include then
+ return next hex;
+ end if;
+
+ -- we could in theory eliminate the inner loop and
+ -- do a return query
+ base = hex.ijk(hex);
+ for r in 1 .. range loop
+ c = base;
+ c[1] := c[1] + r;
+ c[3] := -c[1] - c[2]; -- keep the invariant
+ hex := hex.cantor(c);
+
+ for h in 0 .. r * 6 - 1 loop
+ return next hex;
+ hex = hex.adjacent(h/r+2,hex);
+ end loop;
+ end loop;
+ return;
+end;
+$$ immutable strict;
+
+create function within_array(hex numeric, range integer = 1, include boolean = false) returns numeric[] language 'sql' as $$
+select array_agg(c) from hex.within(hex, range, include) c;
+$$ immutable strict;
+
+create function atrange(hex numeric, range integer = 1) returns setof numeric language 'plpgsql' as $$
+declare
+ c integer[]; -- cantor coordinates
+begin
+ if range = 0 then
+ return next hex;
+ end if;
+ if range <= 0 then
+ return;
+ end if;
+
+ c = ijk(hex);
+ c[1] := c[1] + range;
+ c[3] := -c[1] - c[2]; -- keep the invariant
+ hex := hex.cantor(c);
+
+ for q in 0 .. range * 6 - 1 loop
+ return next hex;
+ hex = hex.adjacent(q/range+2,hex);
+ end loop;
+ return;
+end;
+$$ immutable strict;
+
+select
+xy(ARRAY[1,2]),
+ijk(xy(ARRAY[1,2])),
+xy(ijk(xy(ARRAY[1,2])))
+;
+
+select xy(ARRAY[2,2]), dir, xy(adjhex(dir, ijk(ARRAY[2,2]))) from generate_series(0,5) as dir;
+select
+xy(adjhex(dir, ijk(ARRAY[1,2]))),
+cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))),
+xy(cantor(xy(adjhex(dir, ijk(ARRAY[1,2])))))
+from generate_series(0,5) as dir;
+
+select cantor(ARRAY[4,3]);
+--select * from within(cantor(ARRAY[4,4]), 2);
+select h, xy(h) from within(cantor(ARRAY[4,3]), 2) as h;
projects / hexagon / commitdiff