API type tweaks. Use more doubles internally. Wrap clKernelSetArg.

This is preliminary work for supporting double-precision kernels (issue #13)

CMakeLists.txt

@@ -130,7 +130,7 @@ else()
   # Recommended C++ compilation flags
   # -Weffc++
   set(PROJCL_COMMON_CXX_FLAGS
-	  "-Wextra -Wall -Wno-unused-parameter -Wno-unused-variable -Wpointer-arith -Wcast-align -Wcast-qual -Wfloat-equal -Wredundant-decls -Wno-long-long -fPIC")
+      "-Wextra -Wall -Wno-unused-parameter -Wno-unused-variable -Wpointer-arith -Wcast-align -Wcast-qual -Wfloat-equal -Wredundant-decls -Wno-long-long -fPIC")
 
   if(CMAKE_COMPILER_IS_GNUCC OR CMAKE_COMPILER_IS_GNUCXX)
 

include/projcl/projcl.h

@@ -28,16 +28,16 @@ void pl_unload_code(PLContext *pl_ctx);
 
 
 PLDatumShiftBuffer *pl_load_datum_shift_data(PLContext *pl_ctx, PLSpheroid src_spheroid, 
-                                             const float *xy, int n, int *outError);
+        const float *xy, size_t count, cl_int *outError);
 void pl_unload_datum_shift_data(PLDatumShiftBuffer *pl_buf);
 
 cl_int pl_shift_datum(PLContext *pl_ctx, PLDatum src_datum, PLDatum dst_datum, PLSpheroid dst_spheroid,
-                      PLDatumShiftBuffer *pl_buf, float *xy_out);
+        PLDatumShiftBuffer *pl_buf, float *xy_out);
 
 
 
 
-PLProjectionBuffer *pl_load_projection_data(PLContext *pl_ctx, const float *xy, int n, int copy, int *outError);
+PLProjectionBuffer *pl_load_projection_data(PLContext *pl_ctx, const float *xy, size_t count, cl_bool copy, int *outError);
 void pl_unload_projection_data(PLProjectionBuffer *pl_buf);
 
 
@@ -62,25 +62,25 @@ cl_int pl_project_points_reverse(PLContext *pl_ctx, PLProjection proj, PLProject
 
 
 PLForwardGeodesicFixedDistanceBuffer *pl_load_forward_geodesic_fixed_distance_data(PLContext *pl_ctx,
-    const float *xy_in, int xy_count, const float *az_in, int az_count, cl_int *outError);
+    const float *xy_in, size_t xy_count, const float *az_in, size_t az_count, cl_int *outError);
 void pl_unload_forward_geodesic_fixed_distance_data(PLForwardGeodesicFixedDistanceBuffer *pl_buf);
-cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFixedDistanceBuffer *pl_buf, float *xy_out, 
-    PLSpheroid pl_ell, float distance);
+cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFixedDistanceBuffer *pl_buf,
+        float *xy_out, PLSpheroid pl_ell, double distance);
 
 
 PLForwardGeodesicFixedAngleBuffer *pl_load_forward_geodesic_fixed_angle_data(PLContext *pl_ctx, 
-    const float *dist_in, int dist_count, cl_int *outError);
+    const float *dist_in, size_t dist_count, cl_int *outError);
 void pl_unload_forward_geodesic_fixed_angle_data(PLForwardGeodesicFixedAngleBuffer *pl_buf);
-cl_int pl_forward_geodesic_fixed_angle(PLContext *pl_ctx, PLForwardGeodesicFixedAngleBuffer *pl_buf, float *xy_in, float *xy_out, 
-    PLSpheroid pl_ell, float angle);
+cl_int pl_forward_geodesic_fixed_angle(PLContext *pl_ctx, PLForwardGeodesicFixedAngleBuffer *pl_buf,
+        float *xy_in, float *xy_out, PLSpheroid pl_ell, double angle);
 
 
 
 PLInverseGeodesicBuffer *pl_load_inverse_geodesic_data(PLContext *pl_ctx,
-													   const float *xy1_in, int xy1_count, int xy1_copy,
-													   const float *xy2_in, int xy2_count,
-													   cl_int *outError);
+        const float *xy1_in, size_t xy1_count, cl_bool xy1_copy,
+        const float *xy2_in, size_t xy2_count,
+        cl_int *outError);
 void pl_unload_inverse_geodesic_data(PLInverseGeodesicBuffer *pl_buf);
 
 cl_int pl_inverse_geodesic(PLContext *pl_ctx, PLInverseGeodesicBuffer *pl_buf, float *dist_out,
-						   PLSpheroid pl_ell, float scale);
+        PLSpheroid pl_ell, double scale);

include/projcl/projcl_types.h

@@ -95,18 +95,18 @@ typedef struct PLCode_s {
 
 typedef struct PLSpheroidInfo_s {
 	PLSpheroid   tag;
-	float		major_axis;
-	float		minor_axis;
-	float		ecc;
-	float		ecc2;
-	float		one_ecc2;
-	float		ec;
-	float		inverse_flattening;
-	float		en[8];
-    float       apa[4];
-    float       krueger_A;
-	float		krueger_alpha[8];
-	float		krueger_beta[8];
+	double		major_axis;
+	double		minor_axis;
+	double		ecc;
+	double		ecc2;
+	double		one_ecc2;
+	double		ec;
+	double		inverse_flattening;
+	double		en[8];
+    double       apa[4];
+    double       krueger_A;
+	double		krueger_alpha[8];
+	double		krueger_beta[8];
 } PLSpheroidInfo;
 
 typedef struct PLContext_s {

include/projcl/projcl_warp.h

@@ -84,30 +84,30 @@ typedef struct PLImageArrayBuffer_s {
 PLImageBuffer *pl_load_image(PLContext *pl_ctx,
         cl_channel_order channel_order, cl_channel_type channel_type,
         size_t width, size_t height, size_t row_pitch,
-        const void* pvData, int copy, int *outError);
+        const void* pvData, cl_bool copy, cl_int *outError);
 void pl_unload_image(PLImageBuffer *buf);
 
 PLImageArrayBuffer *pl_load_image_array(PLContext *pl_ctx,
         cl_channel_order channel_order, cl_channel_type channel_type,
         size_t width, size_t height, size_t row_pitch,
         size_t slice_pitch, size_t tiles_across, size_t tiles_down,
-        const void *pvData, int do_copy, cl_int *outError);
+        const void *pvData, cl_bool copy, cl_int *outError);
 void pl_unload_image_array(PLImageArrayBuffer *buf);
 
 cl_int pl_sample_image(PLContext *pl_ctx, PLPointGridBuffer *grid, PLImageBuffer *img,
         PLImageFilter filter, unsigned char *outData);
 cl_int pl_sample_image_array(PLContext *pl_ctx, PLPointGridBuffer *grid, PLImageArrayBuffer *bufs,
         PLImageFilter filter, unsigned char *outData);
 
-PLPointGridBuffer *pl_load_empty_grid(PLContext *pl_ctx, int count_x, int count_y, int *outError);
+PLPointGridBuffer *pl_load_empty_grid(PLContext *pl_ctx, size_t count_x, size_t count_y, int *outError);
 PLPointGridBuffer *pl_load_grid(PLContext *pl_ctx,
-        float origin_x, float width, int count_x,
-        float origin_y, float height, int count_y,
+        double origin_x, double width, size_t count_x,
+        double origin_y, double height, size_t count_y,
         int *outError);
 void pl_unload_grid(PLPointGridBuffer *grid);
 
 cl_int pl_transform_grid(PLContext *pl_ctx, PLPointGridBuffer *src, PLPointGridBuffer *dst,
-        float sx, float sy, float tx, float ty);
+        double sx, double sy, double tx, double ty);
 cl_int pl_shift_grid_datum(PLContext *pl_ctx,
         PLPointGridBuffer *src, PLDatum src_datum, PLSpheroid src_spheroid,
         PLPointGridBuffer *dst, PLDatum dst_datum, PLSpheroid dst_spheroid);

src/CMakeLists.txt

@@ -10,6 +10,8 @@ set(PROJCL_H
 set(PROJCL_C
   projcl_datum.c
   projcl_geodesic.c
+  projcl_kernel.c
+  projcl_kernel.h
   projcl_params.c
   projcl_project.c
   projcl_run.c

src/projcl_datum.c

@@ -9,10 +9,11 @@
 #include <projcl/projcl.h>
 #include <stdlib.h>
 #include "projcl_run.h"
+#include "projcl_kernel.h"
 #include "projcl_util.h"
 
 PLDatumShiftBuffer *pl_load_datum_shift_data(PLContext *pl_ctx, PLSpheroid src_spheroid, 
-                                             const float *xy, int n, cl_int *outError) {
+        const float *xy, size_t n, cl_int *outError) {
     PLDatumShiftBuffer *pl_buf = NULL;
 
     cl_int error = CL_SUCCESS;
@@ -67,7 +68,7 @@ PLDatumShiftBuffer *pl_load_datum_shift_data(PLContext *pl_ctx, PLSpheroid src_s
     pl_buf->y_rw = y_rw;
     pl_buf->z_rw = z_rw;
 
-    error = pl_run_kernel_geodesic_to_cartesian(cartesian_kernel, pl_ctx, pl_buf, src_spheroid);
+    error = pl_run_kernel_geodesic_to_cartesian(pl_ctx, cartesian_kernel, pl_buf, src_spheroid);
 
     pl_buf->xy_out = clCreateBuffer(pl_ctx->ctx, CL_MEM_WRITE_ONLY, 
                                     sizeof(cl_float) * xy_pad_count * 2, NULL, &error);
@@ -119,11 +120,11 @@ cl_int pl_shift_datum(PLContext *pl_ctx, PLDatum src_datum, PLDatum dst_datum, P
     if (geodesic_kernel == NULL || transform_kernel == NULL)
         return CL_INVALID_KERNEL_NAME;
 
-    error = pl_run_kernel_transform_cartesian(transform_kernel, pl_ctx, pl_buf, src_datum, dst_datum);
+    error = pl_run_kernel_transform_cartesian(pl_ctx, transform_kernel, pl_buf, src_datum, dst_datum);
     if (error != CL_SUCCESS)
         return error;
 
-    error = pl_run_kernel_cartesian_to_geodesic(geodesic_kernel, pl_ctx, pl_buf, xy_out, dst_spheroid);
+    error = pl_run_kernel_cartesian_to_geodesic(pl_ctx, geodesic_kernel, pl_buf, xy_out, dst_spheroid);
     if (error != CL_SUCCESS)
         return error;
 

src/projcl_geodesic.c

@@ -10,10 +10,11 @@
 #include <projcl/projcl.h>
 #include "projcl_run.h"
 #include "projcl_util.h"
+#include "projcl_kernel.h"
 #include "projcl_spheroid.h"
 
 PLForwardGeodesicFixedDistanceBuffer *pl_load_forward_geodesic_fixed_distance_data(PLContext *pl_ctx,
-    const float *xy_in, int xy_count, const float *az_in, int az_count, cl_int *outError)
+    const float *xy_in, size_t xy_count, const float *az_in, size_t az_count, cl_int *outError)
 {
 	cl_int error;
 
@@ -106,7 +107,7 @@ void pl_unload_forward_geodesic_fixed_distance_data(PLForwardGeodesicFixedDistan
 }
 
 PLForwardGeodesicFixedAngleBuffer *pl_load_forward_geodesic_fixed_angle_data(PLContext *pl_ctx,
-    const float *dist_in, int dist_count, cl_int *outError) {
+    const float *dist_in, size_t dist_count, cl_int *outError) {
     cl_int error;
 
     PLForwardGeodesicFixedAngleBuffer *pl_buf;
@@ -161,8 +162,8 @@ void pl_unload_forward_geodesic_fixed_angle_data(PLForwardGeodesicFixedAngleBuff
     free(pl_buf);
 }
 
-cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFixedDistanceBuffer *pl_buf, float *xy_out, 
-    PLSpheroid pl_ell, float distance) {
+cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFixedDistanceBuffer *pl_buf,
+        float *xy_out, PLSpheroid pl_ell, double distance) {
 	cl_kernel sincos_kernel = NULL, sincos1_kernel = NULL, fwd_kernel = NULL;
 
 	if (_pl_spheroid_is_spherical(pl_ell)) {
@@ -189,8 +190,8 @@ cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFi
 
 	cl_int error = CL_SUCCESS;
 
-	error |= clSetKernelArg(sincos1_kernel, 0, sizeof(cl_mem), &pl_buf->xy_in);
-	error |= clSetKernelArg(sincos1_kernel, 1, sizeof(cl_mem), &pl_buf->phi_sincos);	
+	error |= pl_set_kernel_arg_mem(pl_ctx, sincos1_kernel, 0, pl_buf->xy_in);
+	error |= pl_set_kernel_arg_mem(pl_ctx, sincos1_kernel, 1, pl_buf->phi_sincos);	
 	if (error != CL_SUCCESS) {
 		return error;
 	}
@@ -201,8 +202,8 @@ cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFi
 		return error;
 	}
 
-	error |= clSetKernelArg(sincos_kernel, 0, sizeof(cl_mem), &pl_buf->az_in);
-	error |= clSetKernelArg(sincos_kernel, 1, sizeof(cl_mem), &pl_buf->az_sincos);
+	error |= pl_set_kernel_arg_mem(pl_ctx, sincos_kernel, 0, pl_buf->az_in);
+	error |= pl_set_kernel_arg_mem(pl_ctx, sincos_kernel, 1, pl_buf->az_sincos);
 	if (error != CL_SUCCESS) {
 		return error;
 	}
@@ -213,13 +214,13 @@ cl_int pl_forward_geodesic_fixed_distance(PLContext *pl_ctx, PLForwardGeodesicFi
 		return error;
 	}
 
-	return pl_run_kernel_forward_geodesic_fixed_distance(fwd_kernel, pl_ctx, pl_buf, 
-                                                         xy_out, pl_ell, distance);
+	return pl_run_kernel_forward_geodesic_fixed_distance(
+            pl_ctx, fwd_kernel, pl_buf, xy_out, pl_ell, distance);
 }
 
-int pl_forward_geodesic_fixed_angle(PLContext *pl_ctx,
-    PLForwardGeodesicFixedAngleBuffer *pl_buf, float *xy_in, float *xy_out,
-    PLSpheroid pl_ell, float angle) {
+int pl_forward_geodesic_fixed_angle(PLContext *pl_ctx, PLForwardGeodesicFixedAngleBuffer *pl_buf,
+        float *xy_in, float *xy_out, PLSpheroid pl_ell, double angle) {
+    double start[2] = { xy_in[0], xy_in[1] };
     cl_kernel fwd_kernel = NULL;
     if (_pl_spheroid_is_spherical(pl_ell)) {
         fwd_kernel = pl_find_kernel(pl_ctx, "pl_forward_geodesic_fixed_angle_s");
@@ -230,14 +231,14 @@ int pl_forward_geodesic_fixed_angle(PLContext *pl_ctx,
         return CL_INVALID_KERNEL_NAME;
     }
 
-    return pl_run_kernel_forward_geodesic_fixed_angle(fwd_kernel, pl_ctx, pl_buf, 
-                                                      xy_in, xy_out, pl_ell, angle);
+    return pl_run_kernel_forward_geodesic_fixed_angle(
+            pl_ctx, fwd_kernel, pl_buf, start, xy_out, pl_ell, angle);
 }
 
 PLInverseGeodesicBuffer *pl_load_inverse_geodesic_data(PLContext *pl_ctx,
-													   const float *xy1_in, int xy1_count, int xy1_copy,
-													   const float *xy2_in, int xy2_count, 
-													   cl_int *outError)
+        const float *xy1_in, size_t xy1_count, cl_bool xy1_copy,
+        const float *xy2_in, size_t xy2_count, 
+        cl_int *outError)
 {
 	cl_int error = CL_SUCCESS;
 
@@ -306,7 +307,7 @@ void pl_unload_inverse_geodesic_data(PLInverseGeodesicBuffer *pl_buf) {
 }
 
 cl_int pl_inverse_geodesic(PLContext *pl_ctx, PLInverseGeodesicBuffer *pl_buf, float *dist_out,
-						   PLSpheroid pl_ell, float scale) {
+        PLSpheroid pl_ell, double scale) {
 	cl_kernel inv_kernel = NULL;
 
 	if (_pl_spheroid_is_spherical(pl_ell)) {
@@ -318,5 +319,5 @@ cl_int pl_inverse_geodesic(PLContext *pl_ctx, PLInverseGeodesicBuffer *pl_buf, f
 		return CL_INVALID_KERNEL_NAME;
 	}
 
-	return pl_run_kernel_inverse_geodesic(inv_kernel, pl_ctx, pl_buf, dist_out, pl_ell, scale);
+	return pl_run_kernel_inverse_geodesic(pl_ctx, inv_kernel, pl_buf, dist_out, pl_ell, scale);
 }

src/projcl_project.c

@@ -14,7 +14,7 @@
 #include <math.h>
 #include <sys/time.h>
 
-PLProjectionBuffer *pl_load_projection_data(PLContext *pl_ctx, const float *xy, int count, int copy, int *outError) {
+PLProjectionBuffer *pl_load_projection_data(PLContext *pl_ctx, const float *xy, size_t count, cl_bool copy, cl_int *outError) {
 	float *xy_pad = NULL;
 	int xy_pad_count = ck_padding(count, PL_FLOAT_VECTOR_SIZE);