Add smoc_disc(order, theta, phi, radius)

We get the function from libhealpix-cxx-dev, while still using
healpix_bare for the basic pixel operations

Author: df7cb

.gitlab-ci.yml

@@ -4,6 +4,8 @@ stages:
 .build: &build
   stage: build
   image: credativ/postgresql-build:${PGVERSION}
+  before_script:
+    - apt-get -y install libhealpix-cxx-dev
   script:
     - make
     - make install

.travis.yml

@@ -19,7 +19,7 @@ dist: xenial
 before_install:
   # extra apt.pg.o.sh options added in version 204, travis currently has 199 (2019-11-27)
   - sudo apt-get -qq update
-  - sudo apt-get -y install postgresql-common
+  - sudo apt-get -y install postgresql-common libhealpix-cxx-dev
 
 install:
   - sudo /usr/share/postgresql-common/pgdg/apt.postgresql.org.sh -p -v $PG_SUPPORTED_VERSIONS -i

Makefile

@@ -63,6 +63,8 @@ else
   include $(top_srcdir)/contrib/contrib-global.mk
 endif
 
+override CPPFLAGS += -I/usr/include/healpix_cxx
+SHLIB_LINK += -lhealpix_cxx
 # link a second time as PGXS does not allow to change the linker
 PGS_LINKER = g++ $(CXXFLAGS) $(filter-out $(CC) $(CFLAGS), $(LINK.shared))
 pgs_link:  $(shlib) $(OBJS) | $(SHLIB_PREREQS)

expected/moc.out

@@ -559,3 +559,33 @@ SELECT intersection(moc) FROM (VALUES ('0/1-4'::smoc), ('0/2-5'), (NULL)) sub(mo
  0/2-4
 (1 row)
 
+SELECT smoc_disc(0, 0, 0, 1);
+ smoc_disc 
+-----------
+ 0/0-3
+(1 row)
+
+SELECT smoc_disc(1, 0, 0, 1);
+      smoc_disc       
+----------------------
+ 1/1-3,5-7,9-11,13-15
+(1 row)
+
+SELECT smoc_disc(2, 0, 0, 1);
+                         smoc_disc                          
+------------------------------------------------------------
+ 1/3,7,11,15 2/5-7,9-11,21-23,25-27,37-39,41-43,53-55,57-59
+(1 row)
+
+SELECT smoc_disc(0, 0, 0, 3.2);
+ smoc_disc 
+-----------
+ 0/0-11
+(1 row)
+
+SELECT smoc_disc(2, 0, 0, 3.2);
+ smoc_disc 
+-----------
+ 0/0-11
+(1 row)
+

moc.c

@@ -22,6 +22,7 @@ PG_FUNCTION_INFO_V1(smoc_superset_spoint);
 PG_FUNCTION_INFO_V1(smoc_not_superset_spoint);
 PG_FUNCTION_INFO_V1(smoc_union);
 PG_FUNCTION_INFO_V1(smoc_intersection);
+PG_FUNCTION_INFO_V1(smoc_disc);
 
 int32 smoc_output_type = 0;
 
@@ -702,3 +703,30 @@ smoc_intersection(PG_FUNCTION_ARGS)
 		PG_RETURN_NULL();
 	}
 }
+
+Datum
+smoc_disc(PG_FUNCTION_ARGS)
+{
+	int		order = PG_GETARG_INT32(0);
+	double	theta = PG_GETARG_FLOAT8(1);
+	double	phi = PG_GETARG_FLOAT8(2);
+	double	radius = PG_GETARG_FLOAT8(3);
+	void*	moc_in_context = create_moc_in_context(moc_error_out);
+	int32	moc_size;
+	Smoc*	moc_ret;
+
+	moc_disc(moc_in_context, order, theta, phi, radius, moc_error_out);
+
+	moc_size = VARHDRSZ + get_moc_size(moc_in_context, moc_error_out);
+	moc_ret = (Smoc*) palloc0(moc_size);
+	SET_VARSIZE(moc_ret, moc_size);
+
+	if (create_moc_release_context(moc_in_context, moc_ret, moc_error_out))
+		PG_RETURN_POINTER(moc_ret);
+	else
+	{
+		ereport(ERROR, (errcode(ERRCODE_INTERNAL_ERROR),
+				errmsg("Internal error in creation of MOC from disc")));
+		PG_RETURN_NULL();
+	}
+}

pgs_moc_ops.sql.in

@@ -63,3 +63,9 @@ CREATE AGGREGATE intersection (smoc) (
 	SFUNC = smoc_intersection,
 	STYPE = smoc
 );
+
+CREATE FUNCTION smoc_disc ("order" int, theta double precision, phi double precision, radius double precision)
+	RETURNS smoc
+	AS 'MODULE_PATHNAME'
+	LANGUAGE C
+	STRICT;

pgs_process_moc.h

@@ -146,6 +146,10 @@ void
 moc_intersection(void* moc_in_context, Smoc* moc_a, int32 moc_a_end, Smoc* moc_b, int32 moc_b_end,
 												pgs_error_handler error_out);
 
+void
+moc_disc(void* moc_in_context, long order, double theta, double phi, double radius,
+												pgs_error_handler error_out);
+
 #ifdef __cplusplus
 }
 #endif

process_moc.cpp

@@ -13,6 +13,8 @@
 #include <exception>
 #include <stdexcept>
 
+#include <healpix_base.h>
+
 #include "pgs_process_moc.h"
 
 #define LAYDEB 0
@@ -921,3 +923,29 @@ moc_intersection(void* moc_in_context, Smoc* moc_a, int32 moc_a_end, Smoc* moc_b
 		}
 	PGS_CATCH
 }
+
+void
+moc_disc(void* moc_in_context, long order, double theta, double phi, double radius,
+												pgs_error_handler error_out)
+{
+	moc_input* p = static_cast<moc_input*>(moc_in_context);
+	moc_input & m = *p;
+	PGS_TRY
+		rangeset<int64> pixset;
+		Healpix_Base2 hp(order, NEST);
+		pointing p(theta, phi);
+
+		hp.query_disc(p, radius, pixset);
+
+		for (tsize j = 0; j < pixset.nranges(); j++)
+		{
+			hpint64 first = pixset.ivbegin(j);
+			hpint64 last = pixset.ivend(j);
+			healpix_convert(first, order); // convert to order 29
+			healpix_convert(last, order);
+			moc_map_entry input(first, last);
+			m.input_map.insert(m.input_map.end(), input);
+		}
+
+	PGS_CATCH
+}

sql/moc.sql

@@ -124,3 +124,9 @@ SELECT smoc_intersection('1/1,4-6', '1/3-5 2/8');
 SELECT '0/1'::smoc * '1/3,5,7,9' AS intersection;
 SELECT '1/9,11,13,15'::smoc * '0/1,2' AS intersection;
 SELECT intersection(moc) FROM (VALUES ('0/1-4'::smoc), ('0/2-5'), (NULL)) sub(moc);
+
+SELECT smoc_disc(0, 0, 0, 1);
+SELECT smoc_disc(1, 0, 0, 1);
+SELECT smoc_disc(2, 0, 0, 1);
+SELECT smoc_disc(0, 0, 0, 3.2);
+SELECT smoc_disc(2, 0, 0, 3.2);