From 7d267065ac29ff006eaf2b316f4a451539d1e55b Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 23 Oct 2024 17:42:29 -0400 Subject: [PATCH] Reformat C code using clang-format --- .clang-format | 17 + src/cdrizzleapi.c | 975 ++++++++++++----------- src/cdrizzleblot.c | 1541 ++++++++++++++++++------------------ src/cdrizzleblot.h | 3 +- src/cdrizzlebox.c | 492 ++++++------ src/cdrizzlebox.h | 14 +- src/cdrizzlemap.h | 91 +-- src/cdrizzleutil.c | 393 +++++---- src/cdrizzleutil.h | 335 ++++---- src/driz_portability.h | 4 +- src/tests/.clang-format | 5 + src/tests/drizzletest.h | 7 +- src/tests/utest_cdrizzle.c | 348 ++++---- 13 files changed, 2095 insertions(+), 2130 deletions(-) create mode 100644 .clang-format create mode 100644 src/tests/.clang-format diff --git a/.clang-format b/.clang-format new file mode 100644 index 00000000..0bee3d2c --- /dev/null +++ b/.clang-format @@ -0,0 +1,17 @@ +--- +BasedOnStyle: Google +Language: Cpp +AlwaysBreakAfterReturnType: None +AlwaysBreakAfterDefinitionReturnType: TopLevel +DerivePointerAlignment: false +PointerAlignment: Right +SortIncludes: false +SortUsingDeclarations: false +IndentWidth: 4 +TabWidth: 4 +UseTab: Never +AlignConsecutiveMacros: true +ReflowComments: true +SpaceBeforeParens: ControlStatements + +... diff --git a/src/cdrizzleapi.c b/src/cdrizzleapi.c index 0fcacb02..4c737088 100644 --- a/src/cdrizzleapi.c +++ b/src/cdrizzleapi.c @@ -1,6 +1,6 @@ #include -#define _USE_MATH_DEFINES /* MS Windows needs to define M_PI */ +#define _USE_MATH_DEFINES /* MS Windows needs to define M_PI */ #include #include #include @@ -23,459 +23,476 @@ static PyObject *gl_Error; FILE *driz_log_handle = NULL; -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Multiply each pixel in an image by a scale factor */ static void scale_image(PyArrayObject *image, double scale_factor) { - long i, size; - float *imptr; + long i, size; + float *imptr; - assert(image); - imptr = (float *) PyArray_DATA(image); + assert(image); + imptr = (float *)PyArray_DATA(image); - size = PyArray_DIMS(image)[0] * PyArray_DIMS(image)[1]; + size = PyArray_DIMS(image)[0] * PyArray_DIMS(image)[1]; - for (i = size; i > 0; --i) { - *imptr++ *= scale_factor; - } + for (i = size; i > 0; --i) { + *imptr++ *= scale_factor; + } - return; + return; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Top level function for drizzling, interfaces with python code */ static PyObject * -tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) -{ - const char *kwlist[] = {"input", "weights", "pixmap", - "output", "counts", "context", - "uniqid", "xmin", "xmax", "ymin", "ymax", - "scale", "pixfrac", "kernel", "in_units", - "expscale", "wtscale", "fillstr", NULL}; - - /* Arguments in the order they appear */ - PyObject *oimg, *owei, *pixmap, *oout, *owht, *ocon; - integer_t uniqid = 1; - integer_t xmin = 0; - integer_t xmax = 0; - integer_t ymin = 0; - integer_t ymax = 0; - double scale = 1.0; - double pfract = 1.0; - char *kernel_str = "square"; - char *inun_str = "cps"; - float expin = 1.0; - float wtscl = 1.0; - char *fillstr = "INDEF"; - - /* Derived values */ - - PyArrayObject *img = NULL, *wei = NULL, *out = NULL, *wht = NULL, *con = NULL, *map = NULL; - enum e_kernel_t kernel; - enum e_unit_t inun; - char *fillstr_end; - bool_t do_fill; - float fill_value; - float inv_exposure_time; - struct driz_error_t error; - struct driz_param_t p; - integer_t isize[2], psize[2], wsize[2]; - char warn_msg[128]; - - driz_log_handle = driz_log_init(driz_log_handle); - driz_log_message("starting tdriz"); - driz_error_init(&error); - - if (!PyArg_ParseTupleAndKeywords(args, keywords, "OOOOOO|iiiiiddssffs:tdriz", (char **)kwlist, - &oimg, &owei, &pixmap, &oout, &owht, &ocon, /* OOOOOO */ - &uniqid, &xmin, &xmax, &ymin, &ymax, /* iiiii */ - &scale, &pfract, &kernel_str, &inun_str, /* ddss */ - &expin, &wtscl, &fillstr) /* ffs */ - ) { - return NULL; - } - - /* Get raw C-array data */ - img = (PyArrayObject *)PyArray_ContiguousFromAny(oimg, NPY_FLOAT, 2, 2); - if (!img) { - driz_error_set_message(&error, "Invalid input array"); - goto _exit; - } - - wei = (PyArrayObject *)PyArray_ContiguousFromAny(owei, NPY_FLOAT, 2, 2); - if (!wei) { - driz_error_set_message(&error, "Invalid weights array"); - goto _exit; - } - - map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); - if (!map) { - driz_error_set_message(&error, "Invalid pixmap array"); - goto _exit; - } - - out = (PyArrayObject *)PyArray_ContiguousFromAny(oout, NPY_FLOAT, 2, 2); - if (!out) { - driz_error_set_message(&error, "Invalid output array"); - goto _exit; - } - - wht = (PyArrayObject *)PyArray_ContiguousFromAny(owht, NPY_FLOAT, 2, 2); - if (!wht) { - driz_error_set_message(&error, "Invalid counts array"); - goto _exit; - } - - if (ocon == Py_None) { - con = NULL; - } else { - con = (PyArrayObject *)PyArray_ContiguousFromAny(ocon, NPY_INT32, 2, 2); - if (!con) { - driz_error_set_message(&error, "Invalid context array"); - goto _exit; - } - }; - - /* Convert the fill value string */ - - if (fillstr == NULL || - *fillstr == 0 || - strncmp(fillstr, "INDEF", 6) == 0 || - strncmp(fillstr, "indef", 6) == 0) { - - do_fill = 0; - fill_value = 0.0; - - } else if (strncmp(fillstr, "NaN", 4) == 0 || - strncmp(fillstr, "nan", 4) == 0) { - do_fill = 1; - fill_value = NPY_NANF; - - } else { - do_fill = 1; +tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords) { + const char *kwlist[] = {"input", "weights", "pixmap", "output", + "counts", "context", "uniqid", "xmin", + "xmax", "ymin", "ymax", "scale", + "pixfrac", "kernel", "in_units", "expscale", + "wtscale", "fillstr", NULL}; + + /* Arguments in the order they appear */ + PyObject *oimg, *owei, *pixmap, *oout, *owht, *ocon; + integer_t uniqid = 1; + integer_t xmin = 0; + integer_t xmax = 0; + integer_t ymin = 0; + integer_t ymax = 0; + double scale = 1.0; + double pfract = 1.0; + char *kernel_str = "square"; + char *inun_str = "cps"; + float expin = 1.0; + float wtscl = 1.0; + char *fillstr = "INDEF"; + + /* Derived values */ + + PyArrayObject *img = NULL, *wei = NULL, *out = NULL, *wht = NULL, + *con = NULL, *map = NULL; + enum e_kernel_t kernel; + enum e_unit_t inun; + char *fillstr_end; + bool_t do_fill; + float fill_value; + float inv_exposure_time; + struct driz_error_t error; + struct driz_param_t p; + integer_t isize[2], psize[2], wsize[2]; + char warn_msg[128]; + + driz_log_handle = driz_log_init(driz_log_handle); + driz_log_message("starting tdriz"); + driz_error_init(&error); + + if (!PyArg_ParseTupleAndKeywords( + args, keywords, "OOOOOO|iiiiiddssffs:tdriz", (char **)kwlist, &oimg, + &owei, &pixmap, &oout, &owht, &ocon, /* OOOOOO */ + &uniqid, &xmin, &xmax, &ymin, &ymax, /* iiiii */ + &scale, &pfract, &kernel_str, &inun_str, /* ddss */ + &expin, &wtscl, &fillstr) /* ffs */ + ) { + return NULL; + } + + /* Get raw C-array data */ + img = (PyArrayObject *)PyArray_ContiguousFromAny(oimg, NPY_FLOAT, 2, 2); + if (!img) { + driz_error_set_message(&error, "Invalid input array"); + goto _exit; + } + + wei = (PyArrayObject *)PyArray_ContiguousFromAny(owei, NPY_FLOAT, 2, 2); + if (!wei) { + driz_error_set_message(&error, "Invalid weights array"); + goto _exit; + } + + map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); + if (!map) { + driz_error_set_message(&error, "Invalid pixmap array"); + goto _exit; + } + + out = (PyArrayObject *)PyArray_ContiguousFromAny(oout, NPY_FLOAT, 2, 2); + if (!out) { + driz_error_set_message(&error, "Invalid output array"); + goto _exit; + } + + wht = (PyArrayObject *)PyArray_ContiguousFromAny(owht, NPY_FLOAT, 2, 2); + if (!wht) { + driz_error_set_message(&error, "Invalid counts array"); + goto _exit; + } + + if (ocon == Py_None) { + con = NULL; + } else { + con = (PyArrayObject *)PyArray_ContiguousFromAny(ocon, NPY_INT32, 2, 2); + if (!con) { + driz_error_set_message(&error, "Invalid context array"); + goto _exit; + } + }; + + /* Convert the fill value string */ + + if (fillstr == NULL || *fillstr == 0 || strncmp(fillstr, "INDEF", 6) == 0 || + strncmp(fillstr, "indef", 6) == 0) { + do_fill = 0; + fill_value = 0.0; + + } else if (strncmp(fillstr, "NaN", 4) == 0 || + strncmp(fillstr, "nan", 4) == 0) { + do_fill = 1; + fill_value = NPY_NANF; + + } else { + do_fill = 1; #ifdef _WIN32 - fill_value = atof(fillstr); + fill_value = atof(fillstr); #else - fill_value = strtof(fillstr, &fillstr_end); - if (fillstr == fillstr_end || *fillstr_end != '\0') { - driz_error_set_message(&error, "Illegal fill value"); - goto _exit; - } + fill_value = strtof(fillstr, &fillstr_end); + if (fillstr == fillstr_end || *fillstr_end != '\0') { + driz_error_set_message(&error, "Illegal fill value"); + goto _exit; + } #endif - } - - /* Set the area to be processed */ - - get_dimensions(img, isize); - if (xmax == 0 || xmax >= isize[0]) xmax = isize[0] - 1; - if (ymax == 0 || ymax >= isize[1]) ymax = isize[1] - 1; - - if (shrink_image_section(map, &xmin, &xmax, &ymin, &ymax)) { - driz_error_set_message(&error, - "No or too few valid pixels in the pixel map."); - goto _exit; - } - - /* Convert strings to enumerations */ - - if (kernel_str2enum(kernel_str, &kernel, &error) || - unit_str2enum(inun_str, &inun, &error)) { - goto _exit; - } - - if (kernel == kernel_gaussian || kernel == kernel_lanczos2 || kernel == kernel_lanczos3) { - if (snprintf(warn_msg, 128, - "Kernel '%s' is not a flux-conserving kernel.", - kernel_str) < 1) { - strcpy(warn_msg, "Selected kernel is not a flux-conserving kernel."); - } - PyErr_WarnEx(PyExc_Warning, warn_msg, 1); - } - - if (pfract <= 0.001){ - printf("kernel reset to POINT due to pfract being set to 0.0...\n"); - kernel_str2enum("point", &kernel, &error); - } - - /* If the input image is not in CPS we need to divide by the exposure */ - if (inun != unit_cps) { - inv_exposure_time = 1.0f / expin; - scale_image(img, inv_exposure_time); - } - - /* Setup reasonable defaults for drizzling */ - driz_param_init(&p); - - p.data = img; - p.weights = wei; - p.pixmap = map; - p.output_data = out; - p.output_counts = wht; - p.output_context = con; - p.uuid = uniqid; - p.xmin = xmin; - p.ymin = ymin; - p.xmax = xmax; - p.ymax = ymax; - p.scale = scale; - p.pixel_fraction = pfract; - p.kernel = kernel; - p.in_units = inun; - p.exposure_time = expin; - p.weight_scale = wtscl; - p.fill_value = fill_value; - p.error = &error; - - if (driz_error_check(&error, "xmin must be >= 0", p.xmin >= 0)) goto _exit; - if (driz_error_check(&error, "ymin must be >= 0", p.ymin >= 0)) goto _exit; - if (driz_error_check(&error, "xmax must be > xmin", p.xmax > p.xmin)) goto _exit; - if (driz_error_check(&error, "ymax must be > ymin", p.ymax > p.ymin)) goto _exit; - if (driz_error_check(&error, "scale must be > 0", p.scale > 0.0)) goto _exit; - if (driz_error_check(&error, "exposure time must be > 0", p.exposure_time)) goto _exit; - if (driz_error_check(&error, "weight scale must be > 0", p.weight_scale > 0.0)) goto _exit; - - get_dimensions(p.pixmap, psize); - if (psize[0] != isize[0] || psize[1] != isize[1]) { - if (snprintf(warn_msg, 128, - "Pixel map dimensions (%d, %d) != input dimensions (%d, %d).", - psize[0], psize[1], isize[0], isize[1]) < 1) { - strcpy(warn_msg, "Pixel map dimensions != input dimensions."); - } - driz_error_set_message(&error, warn_msg); - goto _exit; - } - - if (p.weights) { - get_dimensions(p.weights, wsize); - if (wsize[0] != isize[0] || wsize[1] != isize[1]) { - if (snprintf(warn_msg, 128, - "Weights array dimensions (%d, %d) != input dimensions (%d, %d).", - wsize[0], wsize[1], isize[0], isize[1]) < 1) { - strcpy(warn_msg, "Weights array dimensions != input dimensions."); - } - driz_error_set_message(&error, warn_msg); - goto _exit; - } - } - - if (dobox(&p)) { - goto _exit; - } - - /* Put in the fill values (if defined) */ - if (do_fill) { - put_fill(&p, fill_value); - } - - _exit: - driz_log_message("ending tdriz"); - driz_log_close(driz_log_handle); - Py_XDECREF(con); - Py_XDECREF(img); - Py_XDECREF(wei); - Py_XDECREF(out); - Py_XDECREF(wht); - Py_XDECREF(map); - - if (driz_error_is_set(&error)) { - PyErr_SetString(PyExc_ValueError, driz_error_get_message(&error)); - return NULL; - } else { - return Py_BuildValue("sii", "Callable C-based DRIZZLE Version 1.12 (28th June 2018)", p.nmiss, p.nskip); - } + } + + /* Set the area to be processed */ + + get_dimensions(img, isize); + if (xmax == 0 || xmax >= isize[0]) xmax = isize[0] - 1; + if (ymax == 0 || ymax >= isize[1]) ymax = isize[1] - 1; + + if (shrink_image_section(map, &xmin, &xmax, &ymin, &ymax)) { + driz_error_set_message(&error, + "No or too few valid pixels in the pixel map."); + goto _exit; + } + + /* Convert strings to enumerations */ + + if (kernel_str2enum(kernel_str, &kernel, &error) || + unit_str2enum(inun_str, &inun, &error)) { + goto _exit; + } + + if (kernel == kernel_gaussian || kernel == kernel_lanczos2 || + kernel == kernel_lanczos3) { + if (snprintf(warn_msg, 128, + "Kernel '%s' is not a flux-conserving kernel.", + kernel_str) < 1) { + strcpy(warn_msg, + "Selected kernel is not a flux-conserving kernel."); + } + PyErr_WarnEx(PyExc_Warning, warn_msg, 1); + } + + if (pfract <= 0.001) { + printf("kernel reset to POINT due to pfract being set to 0.0...\n"); + kernel_str2enum("point", &kernel, &error); + } + + /* If the input image is not in CPS we need to divide by the exposure */ + if (inun != unit_cps) { + inv_exposure_time = 1.0f / expin; + scale_image(img, inv_exposure_time); + } + + /* Setup reasonable defaults for drizzling */ + driz_param_init(&p); + + p.data = img; + p.weights = wei; + p.pixmap = map; + p.output_data = out; + p.output_counts = wht; + p.output_context = con; + p.uuid = uniqid; + p.xmin = xmin; + p.ymin = ymin; + p.xmax = xmax; + p.ymax = ymax; + p.scale = scale; + p.pixel_fraction = pfract; + p.kernel = kernel; + p.in_units = inun; + p.exposure_time = expin; + p.weight_scale = wtscl; + p.fill_value = fill_value; + p.error = &error; + + if (driz_error_check(&error, "xmin must be >= 0", p.xmin >= 0)) goto _exit; + if (driz_error_check(&error, "ymin must be >= 0", p.ymin >= 0)) goto _exit; + if (driz_error_check(&error, "xmax must be > xmin", p.xmax > p.xmin)) + goto _exit; + if (driz_error_check(&error, "ymax must be > ymin", p.ymax > p.ymin)) + goto _exit; + if (driz_error_check(&error, "scale must be > 0", p.scale > 0.0)) + goto _exit; + if (driz_error_check(&error, "exposure time must be > 0", p.exposure_time)) + goto _exit; + if (driz_error_check(&error, "weight scale must be > 0", + p.weight_scale > 0.0)) + goto _exit; + + get_dimensions(p.pixmap, psize); + if (psize[0] != isize[0] || psize[1] != isize[1]) { + if (snprintf( + warn_msg, 128, + "Pixel map dimensions (%d, %d) != input dimensions (%d, %d).", + psize[0], psize[1], isize[0], isize[1]) < 1) { + strcpy(warn_msg, "Pixel map dimensions != input dimensions."); + } + driz_error_set_message(&error, warn_msg); + goto _exit; + } + + if (p.weights) { + get_dimensions(p.weights, wsize); + if (wsize[0] != isize[0] || wsize[1] != isize[1]) { + if (snprintf(warn_msg, 128, + "Weights array dimensions (%d, %d) != input " + "dimensions (%d, %d).", + wsize[0], wsize[1], isize[0], isize[1]) < 1) { + strcpy(warn_msg, + "Weights array dimensions != input dimensions."); + } + driz_error_set_message(&error, warn_msg); + goto _exit; + } + } + + if (dobox(&p)) { + goto _exit; + } + + /* Put in the fill values (if defined) */ + if (do_fill) { + put_fill(&p, fill_value); + } + +_exit: + driz_log_message("ending tdriz"); + driz_log_close(driz_log_handle); + Py_XDECREF(con); + Py_XDECREF(img); + Py_XDECREF(wei); + Py_XDECREF(out); + Py_XDECREF(wht); + Py_XDECREF(map); + + if (driz_error_is_set(&error)) { + PyErr_SetString(PyExc_ValueError, driz_error_get_message(&error)); + return NULL; + } else { + return Py_BuildValue( + "sii", "Callable C-based DRIZZLE Version 1.12 (28th June 2018)", + p.nmiss, p.nskip); + } } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Top level function for blotting, interfaces with python code */ static PyObject * -tblot(PyObject *obj, PyObject *args, PyObject *keywords) -{ - const char *kwlist[] = {"source", "pixmap", "output", - "xmin", "xmax", "ymin", "ymax", - "scale", "kscale", "interp", "exptime", - "misval", "sinscl", NULL}; - - /* Arguments in the order they appear */ - PyObject *oimg, *pixmap, *oout; - long xmin = 0; - long xmax = 0; - long ymin = 0; - long ymax = 0; - double scale = 1.0; - float kscale = 1.0; - char *interp_str = "poly5"; - float ef = 1.0; - float misval = 0.0; - float sinscl = 1.0; - - PyArrayObject *img = NULL, *out = NULL, *map = NULL; - enum e_interp_t interp; - int istat = 0; - struct driz_error_t error; - struct driz_param_t p; - integer_t psize[2], osize[2]; - char warn_msg[128]; - - driz_log_handle = driz_log_init(driz_log_handle); - driz_log_message("starting tblot"); - driz_error_init(&error); - - if (!PyArg_ParseTupleAndKeywords(args, keywords, "OOO|lllldfsfff:tblot", (char **)kwlist, - &oimg, &pixmap, &oout, /* OOO */ - &xmin, &xmax, &ymin, &ymax, /* llll */ - &scale, &kscale, &interp_str, &ef, /* dfsf */ - &misval, &sinscl) /* ff */ - ){ - return NULL; - } - - img = (PyArrayObject *)PyArray_ContiguousFromAny(oimg, NPY_FLOAT, 2, 2); - if (!img) { - driz_error_set_message(&error, "Invalid input array"); - goto _exit; - } - - map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); - if (!map) { - driz_error_set_message(&error, "Invalid pixmap array"); - goto _exit; - } - - out = (PyArrayObject *)PyArray_ContiguousFromAny(oout, NPY_FLOAT, 2, 2); - if (!out) { - driz_error_set_message(&error, "Invalid output array"); - goto _exit; - } - - if (interp_str2enum(interp_str, &interp, &error)) { - goto _exit; - } - - get_dimensions(map, psize); - get_dimensions(out, osize); - - if (psize[0] != osize[0] || psize[1] != osize[1]) { - if (snprintf(warn_msg, 128, - "Pixel map dimensions (%d, %d) != output dimensions (%d, %d).", - psize[0], psize[1], osize[0], osize[1]) < 1) { - strcpy(warn_msg, "Pixel map dimensions != output dimensions."); - } - driz_error_set_message(&error, warn_msg); - goto _exit; - } - - if (xmax == 0) xmax = osize[0]; - if (ymax == 0) ymax = osize[1]; - - driz_param_init(&p); - - p.data = img; - p.output_data = out; - p.xmin = xmin; - p.xmax = xmax; - p.ymin = ymin; - p.ymax = ymax; - p.scale = scale; - p.kscale = kscale; - p.in_units = unit_cps; - p.interpolation = interp; - p.ef = ef; - p.misval = misval; - p.sinscl = sinscl; - p.pixmap = map; - p.error = &error; - - if (driz_error_check(&error, "xmin must be >= 0", p.xmin >= 0)) goto _exit; - if (driz_error_check(&error, "ymin must be >= 0", p.ymin >= 0)) goto _exit; - if (driz_error_check(&error, "xmax must be > xmin", p.xmax > p.xmin)) goto _exit; - if (driz_error_check(&error, "ymax must be > ymin", p.ymax > p.ymin)) goto _exit; - if (driz_error_check(&error, "scale must be > 0", p.scale > 0.0)) goto _exit; - if (driz_error_check(&error, "kscale must be > 0", p.kscale > 0.0)) goto _exit; - if (driz_error_check(&error, "exposure time must be > 0", p.ef > 0.0)) goto _exit; - - if (doblot(&p)) goto _exit; - - _exit: - driz_log_message("ending tblot"); - driz_log_close(driz_log_handle); - Py_DECREF(img); - Py_DECREF(out); - Py_DECREF(map); - - if (driz_error_is_set(&error)) { - if (strcmp(driz_error_get_message(&error), "") != 0) - PyErr_SetString(PyExc_Exception, driz_error_get_message(&error)); - return NULL; - } else { - return Py_BuildValue("i",istat); - } -} +tblot(PyObject *obj, PyObject *args, PyObject *keywords) { + const char *kwlist[] = {"source", "pixmap", "output", "xmin", "xmax", + "ymin", "ymax", "scale", "kscale", "interp", + "exptime", "misval", "sinscl", NULL}; + + /* Arguments in the order they appear */ + PyObject *oimg, *pixmap, *oout; + long xmin = 0; + long xmax = 0; + long ymin = 0; + long ymax = 0; + double scale = 1.0; + float kscale = 1.0; + char *interp_str = "poly5"; + float ef = 1.0; + float misval = 0.0; + float sinscl = 1.0; + + PyArrayObject *img = NULL, *out = NULL, *map = NULL; + enum e_interp_t interp; + int istat = 0; + struct driz_error_t error; + struct driz_param_t p; + integer_t psize[2], osize[2]; + char warn_msg[128]; + + driz_log_handle = driz_log_init(driz_log_handle); + driz_log_message("starting tblot"); + driz_error_init(&error); + + if (!PyArg_ParseTupleAndKeywords( + args, keywords, "OOO|lllldfsfff:tblot", (char **)kwlist, &oimg, + &pixmap, &oout, /* OOO */ + &xmin, &xmax, &ymin, &ymax, /* llll */ + &scale, &kscale, &interp_str, &ef, /* dfsf */ + &misval, &sinscl) /* ff */ + ) { + return NULL; + } + + img = (PyArrayObject *)PyArray_ContiguousFromAny(oimg, NPY_FLOAT, 2, 2); + if (!img) { + driz_error_set_message(&error, "Invalid input array"); + goto _exit; + } + + map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); + if (!map) { + driz_error_set_message(&error, "Invalid pixmap array"); + goto _exit; + } + + out = (PyArrayObject *)PyArray_ContiguousFromAny(oout, NPY_FLOAT, 2, 2); + if (!out) { + driz_error_set_message(&error, "Invalid output array"); + goto _exit; + } + if (interp_str2enum(interp_str, &interp, &error)) { + goto _exit; + } + + get_dimensions(map, psize); + get_dimensions(out, osize); -/** -------------------------------------------------------------------------------------------------- + if (psize[0] != osize[0] || psize[1] != osize[1]) { + if (snprintf( + warn_msg, 128, + "Pixel map dimensions (%d, %d) != output dimensions (%d, %d).", + psize[0], psize[1], osize[0], osize[1]) < 1) { + strcpy(warn_msg, "Pixel map dimensions != output dimensions."); + } + driz_error_set_message(&error, warn_msg); + goto _exit; + } + + if (xmax == 0) xmax = osize[0]; + if (ymax == 0) ymax = osize[1]; + + driz_param_init(&p); + + p.data = img; + p.output_data = out; + p.xmin = xmin; + p.xmax = xmax; + p.ymin = ymin; + p.ymax = ymax; + p.scale = scale; + p.kscale = kscale; + p.in_units = unit_cps; + p.interpolation = interp; + p.ef = ef; + p.misval = misval; + p.sinscl = sinscl; + p.pixmap = map; + p.error = &error; + + if (driz_error_check(&error, "xmin must be >= 0", p.xmin >= 0)) goto _exit; + if (driz_error_check(&error, "ymin must be >= 0", p.ymin >= 0)) goto _exit; + if (driz_error_check(&error, "xmax must be > xmin", p.xmax > p.xmin)) + goto _exit; + if (driz_error_check(&error, "ymax must be > ymin", p.ymax > p.ymin)) + goto _exit; + if (driz_error_check(&error, "scale must be > 0", p.scale > 0.0)) + goto _exit; + if (driz_error_check(&error, "kscale must be > 0", p.kscale > 0.0)) + goto _exit; + if (driz_error_check(&error, "exposure time must be > 0", p.ef > 0.0)) + goto _exit; + + if (doblot(&p)) goto _exit; + +_exit: + driz_log_message("ending tblot"); + driz_log_close(driz_log_handle); + Py_DECREF(img); + Py_DECREF(out); + Py_DECREF(map); + + if (driz_error_is_set(&error)) { + if (strcmp(driz_error_get_message(&error), "") != 0) + PyErr_SetString(PyExc_Exception, driz_error_get_message(&error)); + return NULL; + } else { + return Py_BuildValue("i", istat); + } +} + +/** --------------------------------------------------------------------------- * Top level of C unit tests, interfaces with python code */ static PyObject * -test_cdrizzle(PyObject *self, PyObject *args) -{ - PyObject *data, *weights, *pixmap, *output_data, *output_counts, *output_context; - PyArrayObject *dat, *wei, *map, *odat, *ocnt, *ocon; - - int argc = 1; - char *argv[] = {"utest_cdrizzle", NULL}; - - if (!PyArg_ParseTuple(args,"OOOOOO:test_cdrizzle", &data, &weights, &pixmap, - &output_data, &output_counts, &output_context)) { - return NULL; - } - - dat = (PyArrayObject *)PyArray_ContiguousFromAny(data, NPY_FLOAT, 2, 2); - if (! dat) { - return PyErr_Format(gl_Error, "Invalid data array."); - } - - wei = (PyArrayObject *)PyArray_ContiguousFromAny(weights, NPY_FLOAT, 2, 2); - if (! wei) { - return PyErr_Format(gl_Error, "Invalid weghts array."); - } - - map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 2, 4); - if (! map) { - return PyErr_Format(gl_Error, "Invalid pixmap."); - } - - odat = (PyArrayObject *)PyArray_ContiguousFromAny(output_data, NPY_FLOAT, 2, 2); - if (! odat) { - return PyErr_Format(gl_Error, "Invalid output data array."); - } - - ocnt = (PyArrayObject *)PyArray_ContiguousFromAny(output_counts, NPY_FLOAT, 2, 2); - if (! ocnt) { - return PyErr_Format(gl_Error, "Invalid output counts array."); - } - - ocon = (PyArrayObject *)PyArray_ContiguousFromAny(output_context, NPY_INT32, 2, 2); - if (! ocon) { - return PyErr_Format(gl_Error, "Invalid context array"); - } - - set_test_arrays(dat, wei, map, odat, ocnt, ocon); - utest_cdrizzle(argc, argv); - - return Py_BuildValue(""); -} +test_cdrizzle(PyObject *self, PyObject *args) { + PyObject *data, *weights, *pixmap, *output_data, *output_counts, + *output_context; + PyArrayObject *dat, *wei, *map, *odat, *ocnt, *ocon; + + int argc = 1; + char *argv[] = {"utest_cdrizzle", NULL}; + + if (!PyArg_ParseTuple(args, "OOOOOO:test_cdrizzle", &data, &weights, + &pixmap, &output_data, &output_counts, + &output_context)) { + return NULL; + } + dat = (PyArrayObject *)PyArray_ContiguousFromAny(data, NPY_FLOAT, 2, 2); + if (!dat) { + return PyErr_Format(gl_Error, "Invalid data array."); + } + + wei = (PyArrayObject *)PyArray_ContiguousFromAny(weights, NPY_FLOAT, 2, 2); + if (!wei) { + return PyErr_Format(gl_Error, "Invalid weghts array."); + } + + map = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 2, 4); + if (!map) { + return PyErr_Format(gl_Error, "Invalid pixmap."); + } + + odat = (PyArrayObject *)PyArray_ContiguousFromAny(output_data, NPY_FLOAT, 2, + 2); + if (!odat) { + return PyErr_Format(gl_Error, "Invalid output data array."); + } + + ocnt = (PyArrayObject *)PyArray_ContiguousFromAny(output_counts, NPY_FLOAT, + 2, 2); + if (!ocnt) { + return PyErr_Format(gl_Error, "Invalid output counts array."); + } + + ocon = (PyArrayObject *)PyArray_ContiguousFromAny(output_context, NPY_INT32, + 2, 2); + if (!ocon) { + return PyErr_Format(gl_Error, "Invalid context array"); + } + + set_test_arrays(dat, wei, map, odat, ocnt, ocon); + utest_cdrizzle(argc, argv); + + return Py_BuildValue(""); +} static PyObject * -invert_pixmap_wrap(PyObject *self, PyObject *args) -{ +invert_pixmap_wrap(PyObject *self, PyObject *args) { PyObject *pixmap, *xyout, *bbox; PyArrayObject *xyout_arr, *pixmap_arr, *bbox_arr; struct driz_param_t par; @@ -483,20 +500,22 @@ invert_pixmap_wrap(PyObject *self, PyObject *args) npy_intp *ndim, xyin_dim = 2; const double half = 0.5 - DBL_EPSILON; - xyin = (double *) malloc(2 * sizeof(double)); + xyin = (double *)malloc(2 * sizeof(double)); - if (!PyArg_ParseTuple(args,"OOO:invpixmap", &pixmap, &xyout, &bbox)) { - return NULL; + if (!PyArg_ParseTuple(args, "OOO:invpixmap", &pixmap, &xyout, &bbox)) { + return NULL; } - xyout_arr = (PyArrayObject *)PyArray_ContiguousFromAny(xyout, NPY_DOUBLE, 1, 1); + xyout_arr = + (PyArrayObject *)PyArray_ContiguousFromAny(xyout, NPY_DOUBLE, 1, 1); if (!xyout_arr) { - return PyErr_Format(gl_Error, "Invalid xyout array."); + return PyErr_Format(gl_Error, "Invalid xyout array."); } - pixmap_arr = (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); + pixmap_arr = + (PyArrayObject *)PyArray_ContiguousFromAny(pixmap, NPY_DOUBLE, 3, 3); if (!pixmap_arr) { - return PyErr_Format(gl_Error, "Invalid pixmap."); + return PyErr_Format(gl_Error, "Invalid pixmap."); } par.pixmap = pixmap_arr; @@ -508,14 +527,19 @@ invert_pixmap_wrap(PyObject *self, PyObject *args) par.ymin = 0; par.ymax = ndim[0] - 1; } else { - bbox_arr = (PyArrayObject *)PyArray_ContiguousFromAny(bbox, NPY_DOUBLE, 2, 2); + bbox_arr = + (PyArrayObject *)PyArray_ContiguousFromAny(bbox, NPY_DOUBLE, 2, 2); if (!bbox_arr) { - return PyErr_Format(gl_Error, "Invalid input bounding box."); + return PyErr_Format(gl_Error, "Invalid input bounding box."); } - par.xmin = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 0, 0) - half); - par.xmax = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 0, 1) + half); - par.ymin = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 1, 0) - half); - par.ymax = (integer_t)(*(double*) PyArray_GETPTR2(bbox_arr, 1, 1) + half); + par.xmin = + (integer_t)(*(double *)PyArray_GETPTR2(bbox_arr, 0, 0) - half); + par.xmax = + (integer_t)(*(double *)PyArray_GETPTR2(bbox_arr, 0, 1) + half); + par.ymin = + (integer_t)(*(double *)PyArray_GETPTR2(bbox_arr, 1, 0) - half); + par.ymax = + (integer_t)(*(double *)PyArray_GETPTR2(bbox_arr, 1, 1) + half); } xy = (double *)PyArray_DATA(xyout_arr); @@ -525,47 +549,45 @@ invert_pixmap_wrap(PyObject *self, PyObject *args) } PyArrayObject *arr = (PyArrayObject *)PyArray_SimpleNewFromData( - 1, &xyin_dim, NPY_DOUBLE, xyin); + 1, &xyin_dim, NPY_DOUBLE, xyin); PyArray_ENABLEFLAGS(arr, NPY_ARRAY_OWNDATA); return Py_BuildValue("N", arr); } - static PyObject * -clip_polygon_wrap(PyObject *self, PyObject *args) -{ +clip_polygon_wrap(PyObject *self, PyObject *args) { int k; PyObject *pin, *qin; PyArrayObject *pin_arr, *qin_arr; struct polygon p, q, pq; PyObject *list, *tuple; - if (!PyArg_ParseTuple(args,"OO:clip_polygon", &pin, &qin)) { - return NULL; + if (!PyArg_ParseTuple(args, "OO:clip_polygon", &pin, &qin)) { + return NULL; } pin_arr = (PyArrayObject *)PyArray_ContiguousFromAny(pin, NPY_DOUBLE, 2, 2); if (!pin_arr) { - return PyErr_Format(gl_Error, "Invalid P."); + return PyErr_Format(gl_Error, "Invalid P."); } qin_arr = (PyArrayObject *)PyArray_ContiguousFromAny(qin, NPY_DOUBLE, 2, 2); if (!qin_arr) { - return PyErr_Format(gl_Error, "Invalid Q."); + return PyErr_Format(gl_Error, "Invalid Q."); } p.npv = PyArray_SHAPE(pin_arr)[0]; for (k = 0; k < p.npv; ++k) { - p.v[k].x = *((double *) PyArray_GETPTR2(pin_arr, k, 0)); - p.v[k].y = *((double *) PyArray_GETPTR2(pin_arr, k, 1)); + p.v[k].x = *((double *)PyArray_GETPTR2(pin_arr, k, 0)); + p.v[k].y = *((double *)PyArray_GETPTR2(pin_arr, k, 1)); } q.npv = PyArray_SHAPE(qin_arr)[0]; for (k = 0; k < q.npv; ++k) { - q.v[k].x = *((double *) PyArray_GETPTR2(qin_arr, k, 0)); - q.v[k].y = *((double *) PyArray_GETPTR2(qin_arr, k, 1)); + q.v[k].x = *((double *)PyArray_GETPTR2(qin_arr, k, 0)); + q.v[k].y = *((double *)PyArray_GETPTR2(qin_arr, k, 1)); } clip_polygon_to_window(&p, &q, &pq); @@ -582,60 +604,59 @@ clip_polygon_wrap(PyObject *self, PyObject *args) return Py_BuildValue("N", list); } - -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Table of functions callable from python */ static struct PyMethodDef cdrizzle_methods[] = { - {"tdriz", (PyCFunction)tdriz, METH_VARARGS|METH_KEYWORDS, - "tdriz(image, weights, pixmap, output, counts, context, uniqid, xmin, xmax, ymin, ymax, scale, pixfrac, kernel, in_units, expscale, wtscale, fillstr)"}, - {"tblot", (PyCFunction)tblot, METH_VARARGS|METH_KEYWORDS, - "tblot(image, pixmap, output, xmin, xmax, ymin, ymax, scale, kscale, interp, exptime, misval, sinscl)"}, + {"tdriz", (PyCFunction)tdriz, METH_VARARGS | METH_KEYWORDS, + "tdriz(image, weights, pixmap, output, counts, context, uniqid, xmin, " + "xmax, ymin, ymax, scale, pixfrac, kernel, in_units, expscale, wtscale, " + "fillstr)"}, + {"tblot", (PyCFunction)tblot, METH_VARARGS | METH_KEYWORDS, + "tblot(image, pixmap, output, xmin, xmax, ymin, ymax, scale, kscale, " + "interp, exptime, misval, sinscl)"}, {"test_cdrizzle", test_cdrizzle, METH_VARARGS, - "test_cdrizzle(data, weights, pixmap, output_data, output_counts)"}, - {"invert_pixmap", invert_pixmap_wrap, METH_VARARGS, "invert_pixmap(pixmap, xyout, bbox)"}, + "test_cdrizzle(data, weights, pixmap, output_data, output_counts)"}, + {"invert_pixmap", invert_pixmap_wrap, METH_VARARGS, + "invert_pixmap(pixmap, xyout, bbox)"}, {"clip_polygon", clip_polygon_wrap, METH_VARARGS, "clip_polygon(p, q)"}, - {NULL, NULL} /* sentinel */ + {NULL, NULL} /* sentinel */ }; -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- */ #if PY_MAJOR_VERSION < 3 -PyMODINIT_FUNC initcdrizzle(void) -{ +PyMODINIT_FUNC +initcdrizzle(void) { /* Create the module and add the functions */ - (void) Py_InitModule("cdrizzle", cdrizzle_methods); + (void)Py_InitModule("cdrizzle", cdrizzle_methods); /* Check for errors */ - if (PyErr_Occurred()) - Py_FatalError("can't initialize module cdrizzle"); + if (PyErr_Occurred()) Py_FatalError("can't initialize module cdrizzle"); import_array(); } #else -static struct PyModuleDef moduledef = { - PyModuleDef_HEAD_INIT, - "cdrizzle", - NULL, - -1, - cdrizzle_methods, - NULL, - NULL, - NULL, - NULL -}; - -PyMODINIT_FUNC PyInit_cdrizzle(void) -{ +static struct PyModuleDef moduledef = {PyModuleDef_HEAD_INIT, + "cdrizzle", + NULL, + -1, + cdrizzle_methods, + NULL, + NULL, + NULL, + NULL}; + +PyMODINIT_FUNC +PyInit_cdrizzle(void) { PyObject *m; m = PyModule_Create(&moduledef); /* Check for errors */ - if (PyErr_Occurred()) - Py_FatalError("can't initialize module cdrizzle"); + if (PyErr_Occurred()) Py_FatalError("can't initialize module cdrizzle"); import_array(); return m; diff --git a/src/cdrizzleblot.c b/src/cdrizzleblot.c index 5a462cac..56afef01 100644 --- a/src/cdrizzleblot.c +++ b/src/cdrizzleblot.c @@ -7,755 +7,729 @@ #include "cdrizzleutil.h" #include -#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ +#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ #include #include #include #include -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Signature for functions that perform blotting interpolation. */ -typedef int (interp_function)(const void*, - PyArrayObject*, - const float, const float, - /* Output parameters */ - float*, - struct driz_error_t*); +typedef int(interp_function)(const void *, PyArrayObject *, const float, + const float, + /* Output parameters */ + float *, struct driz_error_t *); -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * A standard set of asserts for all of the interpolation functions */ -#define INTERPOLATION_ASSERTS \ - assert(data); \ - assert(isize[0] > 0); \ - assert(isize[1] > 0); \ - assert(x >= 0.0f && x < (float)isize[0]); \ - assert(y >= 0.0f && y < (float)isize[1]); \ - assert(value); \ - assert(error); \ +#define INTERPOLATION_ASSERTS \ + assert(data); \ + assert(isize[0] > 0); \ + assert(isize[1] > 0); \ + assert(x >= 0.0f && x < (float)isize[0]); \ + assert(y >= 0.0f && y < (float)isize[1]); \ + assert(value); \ + assert(error); -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * A structure to hold parameters for sinc interpolation. */ struct sinc_param_t { - /*The scaling factor for sinc interpolation */ - float sinscl; + /*The scaling factor for sinc interpolation */ + float sinscl; }; -/** -------------------------------------------------------------------------------------------------- - * Procedure to evaluate the bicubic polynomial interpolant. The array coeff contains the coefficients - * of the 2D interpolant. The procedure assumes that 1 <= x <= isize[0] and 1 <= y <= isize[1] and - * that coeff[1+first_point] = datain[1,1]. The interpolant is evaluated using Everett's central - * difference formula. (Was: IIBIP3) +/** --------------------------------------------------------------------------- + * Procedure to evaluate the bicubic polynomial interpolant. The array coeff + * contains the coefficients of the 2D interpolant. The procedure assumes that + * 1 <= x <= isize[0] and 1 <= y <= isize[1] and that coeff[1+first_point] = + * datain[1,1]. The interpolant is evaluated using Everett's central difference + * formula. (Was: IIBIP3) * - * coeff: Array of shape \a [len_coeff][len_coeff] contains the coefficients of the 2D interpolant. - * len_coeff: The dimension (each side of the square) of the coefficient array. - * firstt: Offset of the first data point. (In practice, this is always zero.) - * npts: The number of points to calculate. + * coeff: Array of shape \a [len_coeff][len_coeff] contains the coefficients + * of the 2D interpolant. len_coeff: The dimension (each side of the square) of + * the coefficient array. firstt: Offset of the first data point. (In + * practice, this is always zero.) npts: The number of points to calculate. * x: An array of length \a npts of x values. * y: An array of length \a npts of y values. * zfit: An array of length \a npts of interpolated values. (output) */ static inline_macro void -ii_bipoly3(const float* coeff /* [len_coeff][len_coeff] */, +ii_bipoly3(const float *coeff /* [len_coeff][len_coeff] */, const integer_t len_coeff, const integer_t firstt, - const integer_t npts, - const float* x /* [npts] */, const float* y /* [npts] */, + const integer_t npts, const float *x /* [npts] */, + const float *y /* [npts] */, /* Output parameters */ - float* zfit /* [npts] */) { - float sx, tx, sx2m1, tx2m1, sy, ty; - float cd20[4], cd21[4], ztemp[4]; - float cd20y, cd21y; - integer_t nxold, nyold; - integer_t nx, ny; - integer_t firstw, index; - integer_t i, j; - - nxold = nyold = -1; - for (i = 0; i < npts; ++i) { - nx = (integer_t)x[i]; - assert(nx >= 0); - - sx = x[i] - (float)nx; - tx = 1.0f - sx; - sx2m1 = sx*sx - 1.0f; - tx2m1 = tx*tx - 1.0f; - - ny = (integer_t)y[i]; - assert(ny >= 0); - - sy = y[i] - (float)ny; - ty = 1.0f - sy; - - /* Calculate pointer to data[nx, ny-1] */ - firstw = firstt + (ny - 2) * len_coeff + nx - 1; - - /* loop over the 4 surrounding rows of data calculate the central - differences at each value of y - - If new data point calculate the central differences in x for - each y */ - if (nx != nxold || ny != nyold) { - for (j = 0, index = firstw; j < 4; ++j, index += len_coeff) { - assert(index > 0 && index < (len_coeff*len_coeff) - 2); - - cd20[j] = 1.0f/6.0f * (coeff[index+1] - - 2.0f * coeff[index] + - coeff[index-1]); - cd21[j] = 1.0f/6.0f * (coeff[index+2] - - 2.0f * coeff[index+1] + - coeff[index]); - } - } + float *zfit /* [npts] */) { + float sx, tx, sx2m1, tx2m1, sy, ty; + float cd20[4], cd21[4], ztemp[4]; + float cd20y, cd21y; + integer_t nxold, nyold; + integer_t nx, ny; + integer_t firstw, index; + integer_t i, j; + + nxold = nyold = -1; + for (i = 0; i < npts; ++i) { + nx = (integer_t)x[i]; + assert(nx >= 0); + + sx = x[i] - (float)nx; + tx = 1.0f - sx; + sx2m1 = sx * sx - 1.0f; + tx2m1 = tx * tx - 1.0f; + + ny = (integer_t)y[i]; + assert(ny >= 0); + + sy = y[i] - (float)ny; + ty = 1.0f - sy; + + /* Calculate pointer to data[nx, ny-1] */ + firstw = firstt + (ny - 2) * len_coeff + nx - 1; + + /* loop over the 4 surrounding rows of data calculate the central + differences at each value of y + + If new data point calculate the central differences in x for + each y */ + if (nx != nxold || ny != nyold) { + for (j = 0, index = firstw; j < 4; ++j, index += len_coeff) { + assert(index > 0 && index < (len_coeff * len_coeff) - 2); + + cd20[j] = + 1.0f / 6.0f * + (coeff[index + 1] - 2.0f * coeff[index] + coeff[index - 1]); + cd21[j] = + 1.0f / 6.0f * + (coeff[index + 2] - 2.0f * coeff[index + 1] + coeff[index]); + } + } - /* Interpolate in x at each value of y */ - for (j = 0, index = firstw; j < 4; ++j, index += len_coeff) { - assert(index >= 0 && index < (len_coeff*len_coeff) - 1); + /* Interpolate in x at each value of y */ + for (j = 0, index = firstw; j < 4; ++j, index += len_coeff) { + assert(index >= 0 && index < (len_coeff * len_coeff) - 1); - ztemp[j] = sx * (coeff[index+1] + sx2m1 * cd21[j]) + - tx * (coeff[index] + tx2m1 * cd20[j]); - } + ztemp[j] = sx * (coeff[index + 1] + sx2m1 * cd21[j]) + + tx * (coeff[index] + tx2m1 * cd20[j]); + } - /* Calculate y central differences */ - cd20y = 1.0f/6.0f * (ztemp[2] - 2.0f * ztemp[1] + ztemp[0]); - cd21y = 1.0f/6.0f * (ztemp[3] - 2.0f * ztemp[2] + ztemp[1]); + /* Calculate y central differences */ + cd20y = 1.0f / 6.0f * (ztemp[2] - 2.0f * ztemp[1] + ztemp[0]); + cd21y = 1.0f / 6.0f * (ztemp[3] - 2.0f * ztemp[2] + ztemp[1]); - /* Interpolate in y */ - zfit[i] = sy * (ztemp[2] + (sy * sy - 1.0f) * cd21y) + - ty * (ztemp[1] + (ty * ty - 1.0f) * cd20y); + /* Interpolate in y */ + zfit[i] = sy * (ztemp[2] + (sy * sy - 1.0f) * cd21y) + + ty * (ztemp[1] + (ty * ty - 1.0f) * cd20y); - nxold = nx; - nyold = ny; - } + nxold = nx; + nyold = ny; + } } -/** -------------------------------------------------------------------------------------------------- - * Procedure to evaluate a biquintic polynomial. The array coeff contains the coefficents of the - * 2D interpolant. The routine assumes that 0 <= x < isize[0] and 0 <= y < isize[1]. The interpolant - * is evaluated using Everett's central difference formula. (Was: IIBIP5) +/** --------------------------------------------------------------------------- + * Procedure to evaluate a biquintic polynomial. The array coeff contains the + * coefficents of the 2D interpolant. The routine assumes that 0 <= x < + * isize[0] and 0 <= y < isize[1]. The interpolant is evaluated using Everett's + * central difference formula. (Was: IIBIP5) * - * coeff: Array of shape \a [len_coeff][len_coeff] contains the coefficients of the 2D interpolant. - * len_coeff: The dimension (each side of the square) of the coefficient array. - * firstt: Offset of the first data point. (In practice, this is always zero.) - * npts: The number of points to calculate. + * coeff: Array of shape \a [len_coeff][len_coeff] contains the coefficients + * of the 2D interpolant. len_coeff: The dimension (each side of the square) of + * the coefficient array. firstt: Offset of the first data point. (In + * practice, this is always zero.) npts: The number of points to calculate. * x: An array of length \a npts of x values. * y: An array of length \a npts of y values. * zfit: An array of length \a npts of interpolated values. (output) */ static inline_macro void -ii_bipoly5(const float* coeff /* [len_coeff][len_coeff] */, +ii_bipoly5(const float *coeff /* [len_coeff][len_coeff] */, const integer_t len_coeff, const integer_t firstt, - const integer_t npts, - const float* x /* [npts] */, const float* y /* [npts] */, + const integer_t npts, const float *x /* [npts] */, + const float *y /* [npts] */, /* Output parameters */ - float* zfit /* [npts] */) { - integer_t nxold, nyold; - integer_t nx, ny; - float sx, sx2, tx, tx2, sy, sy2, ty, ty2; - float sx2m1, sx2m4, tx2m1, tx2m4; - float cd20[6], cd21[6], cd40[6], cd41[6]; - float cd20y, cd21y, cd40y, cd41y; - float ztemp[6]; - integer_t firstw, index; - integer_t i, j; - - assert(coeff); - assert(len_coeff > 0); - assert(npts > 0); - assert(x); - assert(y); - assert(zfit); - - nxold = nyold = -1; - for (i = 0; i < npts; ++i) { - nx = (integer_t)x[i]; - ny = (integer_t)y[i]; - assert(nx >= 0); - assert(ny >= 0); - - sx = x[i] - (float)nx; - sx2 = sx * sx; - sx2m1 = sx2 - 1.0f; - sx2m4 = sx2 - 4.0f; - tx = 1.0f - sx; - tx2 = tx * tx; - tx2m1 = tx2 - 1.0f; - tx2m4 = tx2 - 4.0f; - - sy = y[i] - (float)ny; - sy2 = sy * sy; - ty = 1.0f - sy; - ty2 = ty * ty; - - /* Calculate value of pointer to data[nx,ny-2] */ - firstw = firstt + (ny - 3)*len_coeff + nx - 1; - - /* Calculate the central differences in x at each value of y */ - if (nx != nxold || ny != nyold) { - for (j = 0, index = firstw; j < 6; ++j, index += len_coeff) { - assert(index >= 1 && index < len_coeff*len_coeff - 3); - - cd20[j] = 1.0f/6.0f * (coeff[index+1] - - 2.0f * coeff[index] + - coeff[index-1]); - cd21[j] = 1.0f/6.0f * (coeff[index+2] - - 2.0f * coeff[index+1] + - coeff[index]); - cd40[j] = 1.0f/120.0f * (coeff[index-2] - - 4.0f * coeff[index-1] + - 6.0f * coeff[index] - - 4.0f * coeff[index+1] + - coeff[index+2]); - cd41[j] = 1.0f/120.0f * (coeff[index-1] - - 4.0f * coeff[index] + - 6.0f * coeff[index+1] - - 4.0f * coeff[index+2] + - coeff[index+3]); - } - } + float *zfit /* [npts] */) { + integer_t nxold, nyold; + integer_t nx, ny; + float sx, sx2, tx, tx2, sy, sy2, ty, ty2; + float sx2m1, sx2m4, tx2m1, tx2m4; + float cd20[6], cd21[6], cd40[6], cd41[6]; + float cd20y, cd21y, cd40y, cd41y; + float ztemp[6]; + integer_t firstw, index; + integer_t i, j; + + assert(coeff); + assert(len_coeff > 0); + assert(npts > 0); + assert(x); + assert(y); + assert(zfit); + + nxold = nyold = -1; + for (i = 0; i < npts; ++i) { + nx = (integer_t)x[i]; + ny = (integer_t)y[i]; + assert(nx >= 0); + assert(ny >= 0); + + sx = x[i] - (float)nx; + sx2 = sx * sx; + sx2m1 = sx2 - 1.0f; + sx2m4 = sx2 - 4.0f; + tx = 1.0f - sx; + tx2 = tx * tx; + tx2m1 = tx2 - 1.0f; + tx2m4 = tx2 - 4.0f; + + sy = y[i] - (float)ny; + sy2 = sy * sy; + ty = 1.0f - sy; + ty2 = ty * ty; + + /* Calculate value of pointer to data[nx,ny-2] */ + firstw = firstt + (ny - 3) * len_coeff + nx - 1; + + /* Calculate the central differences in x at each value of y */ + if (nx != nxold || ny != nyold) { + for (j = 0, index = firstw; j < 6; ++j, index += len_coeff) { + assert(index >= 1 && index < len_coeff * len_coeff - 3); + + cd20[j] = + 1.0f / 6.0f * + (coeff[index + 1] - 2.0f * coeff[index] + coeff[index - 1]); + cd21[j] = + 1.0f / 6.0f * + (coeff[index + 2] - 2.0f * coeff[index + 1] + coeff[index]); + cd40[j] = 1.0f / 120.0f * + (coeff[index - 2] - 4.0f * coeff[index - 1] + + 6.0f * coeff[index] - 4.0f * coeff[index + 1] + + coeff[index + 2]); + cd41[j] = 1.0f / 120.0f * + (coeff[index - 1] - 4.0f * coeff[index] + + 6.0f * coeff[index + 1] - 4.0f * coeff[index + 2] + + coeff[index + 3]); + } + } - /* Interpolate in x at each value of y */ - for (j = 0, index = firstw; j < 6; ++j, index += len_coeff) { - assert(index >= 0 && index < len_coeff*len_coeff - 1); + /* Interpolate in x at each value of y */ + for (j = 0, index = firstw; j < 6; ++j, index += len_coeff) { + assert(index >= 0 && index < len_coeff * len_coeff - 1); - ztemp[j] = sx * (coeff[index+1] + sx2m1 * (cd21[j] + sx2m4 * cd41[j])) + - tx * (coeff[index] + tx2m1 * (cd20[j] + tx2m4 * cd40[j])); - } + ztemp[j] = + sx * (coeff[index + 1] + sx2m1 * (cd21[j] + sx2m4 * cd41[j])) + + tx * (coeff[index] + tx2m1 * (cd20[j] + tx2m4 * cd40[j])); + } - /* Central differences in y */ - cd20y = 1.0f/6.0f * (ztemp[3] - 2.0f * ztemp[2] + ztemp[1]); - cd21y = 1.0f/6.0f * (ztemp[4] - 2.0f * ztemp[3] + ztemp[2]); - cd40y = 1.0f/120.0f * (ztemp[0] - - 4.0f * ztemp[1] + - 6.0f * ztemp[2] - - 4.0f * ztemp[3] + - ztemp[4]); - cd41y = 1.0f/120.0f * (ztemp[1] - - 4.0f * ztemp[2] + - 6.0f * ztemp[3] - - 4.0f * ztemp[4] + - ztemp[5]); - - /* Interpolate in y */ - zfit[i] = sy * (ztemp[3] + (sy2 - 1.0f) * (cd21y + (sy2 - 4.0f) * cd41y)) + - ty * (ztemp[2] + (ty2 - 1.0f) * (cd20y + (ty2 - 4.0f) * cd40y)); - - nxold = nx; - nyold = ny; - } + /* Central differences in y */ + cd20y = 1.0f / 6.0f * (ztemp[3] - 2.0f * ztemp[2] + ztemp[1]); + cd21y = 1.0f / 6.0f * (ztemp[4] - 2.0f * ztemp[3] + ztemp[2]); + cd40y = 1.0f / 120.0f * + (ztemp[0] - 4.0f * ztemp[1] + 6.0f * ztemp[2] - + 4.0f * ztemp[3] + ztemp[4]); + cd41y = 1.0f / 120.0f * + (ztemp[1] - 4.0f * ztemp[2] + 6.0f * ztemp[3] - + 4.0f * ztemp[4] + ztemp[5]); + + /* Interpolate in y */ + zfit[i] = + sy * (ztemp[3] + (sy2 - 1.0f) * (cd21y + (sy2 - 4.0f) * cd41y)) + + ty * (ztemp[2] + (ty2 - 1.0f) * (cd20y + (ty2 - 4.0f) * cd40y)); + + nxold = nx; + nyold = ny; + } } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform nearest neighbor interpolation. * - * state: A pointer to any constant values specific to this interpolation type. (NULL). - * data: A 2D data array - * x: The fractional x coordinate - * y: The fractional y coordinate - * value: The resulting value at x, y after interpolating the data (output) - * error: The error structure (output) + * state: A pointer to any constant values specific to this interpolation type. + * (NULL). data: A 2D data array x: The fractional x coordinate y: The + * fractional y coordinate value: The resulting value at x, y after + * interpolating the data (output) error: The error structure (output) */ static int -interpolate_nearest_neighbor(const void* state UNUSED_PARAM, - PyArrayObject* data, - const float x, const float y, +interpolate_nearest_neighbor(const void *state UNUSED_PARAM, + PyArrayObject *data, const float x, const float y, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - - integer_t isize[2]; - get_dimensions(data, isize); + float *value, + struct driz_error_t *error UNUSED_PARAM) { + integer_t isize[2]; + get_dimensions(data, isize); - assert(state == NULL); - INTERPOLATION_ASSERTS; + assert(state == NULL); + INTERPOLATION_ASSERTS; - *value = get_pixel(data, (integer_t)(x + 0.5), (integer_t)(y + 0.5)); - return 0; + *value = get_pixel(data, (integer_t)(x + 0.5), (integer_t)(y + 0.5)); + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform basic bilinear interpolation. * - * state: A pointer to any constant values specific to this interpolation type. (NULL). - * data: A 2D data array - * x: The fractional x coordinate - * y: The fractional y coordinate - * value: The resulting value at x, y after interpolating the data (output) - * error: The error structure (output) + * state: A pointer to any constant values specific to this interpolation type. + * (NULL). data: A 2D data array x: The fractional x coordinate y: The + * fractional y coordinate value: The resulting value at x, y after + * interpolating the data (output) error: The error structure (output) */ static int -interpolate_bilinear(const void* state UNUSED_PARAM, - PyArrayObject* data, +interpolate_bilinear(const void *state UNUSED_PARAM, PyArrayObject *data, const float x, const float y, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - integer_t nx, ny; - float sx, tx, sy, ty, f00; - integer_t isize[2]; + float *value, struct driz_error_t *error UNUSED_PARAM) { + integer_t nx, ny; + float sx, tx, sy, ty, f00; + integer_t isize[2]; - get_dimensions(data, isize); + get_dimensions(data, isize); - assert(state == NULL); - INTERPOLATION_ASSERTS; + assert(state == NULL); + INTERPOLATION_ASSERTS; - nx = (integer_t) x; - ny = (integer_t) y; + nx = (integer_t)x; + ny = (integer_t)y; - if (nx < 0 || ny < 0 || nx >= isize[0] || ny >= isize[1]) { - driz_error_set_message(error, - "Bilinear interpolation: point outside of the image."); - return 1; - } + if (nx < 0 || ny < 0 || nx >= isize[0] || ny >= isize[1]) { + driz_error_set_message( + error, "Bilinear interpolation: point outside of the image."); + return 1; + } - f00 = get_pixel(data, nx, ny); + f00 = get_pixel(data, nx, ny); - if (nx == (isize[0] - 1)) { - if (ny == (isize[1] - 1)) { - /* This is the last pixel (in x and y). Assign constant value of this pixel. */ - *value = f00; - return 0; + if (nx == (isize[0] - 1)) { + if (ny == (isize[1] - 1)) { + /* This is the last pixel (in x and y). Assign constant value of + * this pixel. */ + *value = f00; + return 0; + } + /* Interpolate along Y-direction only */ + sy = y - (float)ny; + *value = (1.0f - sy) * f00 + sy * get_pixel(data, nx, ny + 1); + } else if (ny == (isize[1] - 1)) { + /* Interpolate along X-direction only */ + sx = x - (float)nx; + *value = (1.0f - sx) * f00 + sx * get_pixel(data, nx + 1, ny); + } else { + /* Bilinear - interpolation */ + sx = x - (float)nx; + tx = 1.0f - sx; + sy = y - (float)ny; + ty = 1.0f - sy; + + *value = tx * ty * f00 + sx * ty * get_pixel(data, nx + 1, ny) + + sy * tx * get_pixel(data, nx, ny + 1) + + sx * sy * get_pixel(data, nx + 1, ny + 1); } - /* Interpolate along Y-direction only */ - sy = y - (float)ny; - *value = (1.0f - sy) * f00 + sy * get_pixel(data, nx, ny + 1); - } else if (ny == (isize[1] - 1)) { - /* Interpolate along X-direction only */ - sx = x - (float)nx; - *value = (1.0f - sx) * f00 + sx * get_pixel(data, nx + 1, ny); - } else { - /* Bilinear - interpolation */ - sx = x - (float)nx; - tx = 1.0f - sx; - sy = y - (float)ny; - ty = 1.0f - sy; - - *value = tx * ty * f00 + - sx * ty * get_pixel(data, nx + 1, ny) + - sy * tx * get_pixel(data, nx, ny + 1) + - sx * sy * get_pixel(data, nx + 1, ny + 1); - } - - return 0; + + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform cubic polynomial interpolation. * - * state: A pointer to any constant values specific to this interpolation type. (NULL). - * data: A 2D data array - * x: The fractional x coordinate - * y: The fractional y coordinate - * value: The resulting value at x, y after interpolating the data (output) - * error: The error structure (output) + * state: A pointer to any constant values specific to this interpolation type. + * (NULL). data: A 2D data array x: The fractional x coordinate y: The + * fractional y coordinate value: The resulting value at x, y after + * interpolating the data (output) error: The error structure (output) */ static int -interpolate_poly3(const void* state UNUSED_PARAM, - PyArrayObject* data, +interpolate_poly3(const void *state UNUSED_PARAM, PyArrayObject *data, const float x, const float y, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - integer_t nx, ny; - const integer_t rowleh = 4; - const integer_t nterms = 4; - float coeff[4][4]; - integer_t i, j; - integer_t firstw, lastrw; - float xval, yval; - float* ci; - integer_t isize[2]; - get_dimensions(data, isize); - - assert(state == NULL); - INTERPOLATION_ASSERTS;; - - nx = (integer_t)x; - ny = (integer_t)y; - - ci = &coeff[0][0]; - for (j = ny - 1; j <= ny + 2; ++j) { - if (j >= 0 && j < isize[1]) { - for (i = nx - 1; i <= nx + 2; ++i, ++ci) { - if (i < 0) { - *ci = 2.0f * get_pixel(data, 0, j) - get_pixel(data, -i, j); - } else if (i >= isize[0]) { - *ci = 2.0f * get_pixel(data, isize[0]-1, j) - get_pixel(data, 2*isize[0]-2-i, j); - } else { - *ci = get_pixel(data, i, j); - } - } - } else if (j == ny + 2) { - for (i = nx - 1; i <= nx + 2; ++i, ++ci) { - if (i < 0) { - *ci = 2.0f * get_pixel(data, 0, isize[1]-3) - get_pixel(data, -i, isize[1]-3); - } else if (i >= isize[0]) { - *ci = 2.0f * get_pixel(data, isize[0]-1, isize[1]-3) - get_pixel(data, 2*isize[0]-2-i, isize[1]-3); + float *value, struct driz_error_t *error UNUSED_PARAM) { + integer_t nx, ny; + const integer_t rowleh = 4; + const integer_t nterms = 4; + float coeff[4][4]; + integer_t i, j; + integer_t firstw, lastrw; + float xval, yval; + float *ci; + integer_t isize[2]; + get_dimensions(data, isize); + + assert(state == NULL); + INTERPOLATION_ASSERTS; + ; + + nx = (integer_t)x; + ny = (integer_t)y; + + ci = &coeff[0][0]; + for (j = ny - 1; j <= ny + 2; ++j) { + if (j >= 0 && j < isize[1]) { + for (i = nx - 1; i <= nx + 2; ++i, ++ci) { + if (i < 0) { + *ci = 2.0f * get_pixel(data, 0, j) - get_pixel(data, -i, j); + } else if (i >= isize[0]) { + *ci = 2.0f * get_pixel(data, isize[0] - 1, j) - + get_pixel(data, 2 * isize[0] - 2 - i, j); + } else { + *ci = get_pixel(data, i, j); + } + } + } else if (j == ny + 2) { + for (i = nx - 1; i <= nx + 2; ++i, ++ci) { + if (i < 0) { + *ci = 2.0f * get_pixel(data, 0, isize[1] - 3) - + get_pixel(data, -i, isize[1] - 3); + } else if (i >= isize[0]) { + *ci = 2.0f * get_pixel(data, isize[0] - 1, isize[1] - 3) - + get_pixel(data, 2 * isize[0] - 2 - i, isize[1] - 3); + } else { + *ci = get_pixel(data, i, isize[1] - 3); + } + } } else { - *ci = get_pixel(data, i, isize[1]-3); + ci += 4; } - } - } else { - ci += 4; } - } - firstw = MAX(0, 1 - ny); - if (firstw > 0) { - assert(firstw < nterms); + firstw = MAX(0, 1 - ny); + if (firstw > 0) { + assert(firstw < nterms); - for (j = 0; j < firstw; ++j) { - assert(2*firstw-j >= 0 && 2*firstw-j < nterms); + for (j = 0; j < firstw; ++j) { + assert(2 * firstw - j >= 0 && 2 * firstw - j < nterms); - weighted_sum_vectors(nterms, - &coeff[firstw][0], 2.0, - &coeff[2*firstw-j][0], -1.0, - &coeff[j][0]); + weighted_sum_vectors(nterms, &coeff[firstw][0], 2.0, + &coeff[2 * firstw - j][0], -1.0, &coeff[j][0]); + } } - } - lastrw = MIN(nterms - 1, isize[1] - ny); - if (lastrw < nterms - 1) { - assert(lastrw >= 0 && lastrw < nterms); + lastrw = MIN(nterms - 1, isize[1] - ny); + if (lastrw < nterms - 1) { + assert(lastrw >= 0 && lastrw < nterms); - for (j = lastrw + 1; j <= nterms - 1; ++j) { - assert(2*lastrw-j >= 0 && 2*lastrw-j < nterms); - assert(j >= 0 && j < 4); + for (j = lastrw + 1; j <= nterms - 1; ++j) { + assert(2 * lastrw - j >= 0 && 2 * lastrw - j < nterms); + assert(j >= 0 && j < 4); - weighted_sum_vectors(nterms, - &coeff[lastrw][0], 2.0, - &coeff[2*lastrw-j][0], -1.0, - &coeff[j][0]); - } - } else if (lastrw == 1) { - assert(lastrw >= 0 && lastrw < nterms); + weighted_sum_vectors(nterms, &coeff[lastrw][0], 2.0, + &coeff[2 * lastrw - j][0], -1.0, &coeff[j][0]); + } + } else if (lastrw == 1) { + assert(lastrw >= 0 && lastrw < nterms); - weighted_sum_vectors(nterms, - &coeff[lastrw][0], 2.0, - &coeff[3][0], -1.0, - &coeff[3][0]); - } else { - assert(lastrw >= 0 && lastrw < nterms); - assert(2*lastrw-3 >= 0 && 2*lastrw-3 < nterms); + weighted_sum_vectors(nterms, &coeff[lastrw][0], 2.0, &coeff[3][0], -1.0, + &coeff[3][0]); + } else { + assert(lastrw >= 0 && lastrw < nterms); + assert(2 * lastrw - 3 >= 0 && 2 * lastrw - 3 < nterms); - weighted_sum_vectors(nterms, - &coeff[lastrw][0], 2.0, - &coeff[2*lastrw-3][0], -1.0, - &coeff[3][0]); - } + weighted_sum_vectors(nterms, &coeff[lastrw][0], 2.0, + &coeff[2 * lastrw - 3][0], -1.0, &coeff[3][0]); + } - xval = 2.0f + (x - (float)nx); - yval = 2.0f + (y - (float)ny); + xval = 2.0f + (x - (float)nx); + yval = 2.0f + (y - (float)ny); - ii_bipoly3(&coeff[0][0], rowleh, 0, 1, &xval, &yval, value); + ii_bipoly3(&coeff[0][0], rowleh, 0, 1, &xval, &yval, value); - return 0; + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform quintic polynomial interpolation. * - * state: A pointer to any constant values specific to this interpolation type. (NULL). - * data: A 2D data array - * x: The fractional x coordinate - * y: The fractional y coordinate - * value: The resulting value at x, y after interpolating the data (output) - * error: The error structure (output) + * state: A pointer to any constant values specific to this interpolation type. + * (NULL). data: A 2D data array x: The fractional x coordinate y: The + * fractional y coordinate value: The resulting value at x, y after + * interpolating the data (output) error: The error structure (output) */ static int -interpolate_poly5(const void* state UNUSED_PARAM, - PyArrayObject* data, +interpolate_poly5(const void *state UNUSED_PARAM, PyArrayObject *data, const float x, const float y, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - integer_t nx, ny; - const integer_t rowleh = 6; - const integer_t nterms = 6; - float coeff[6][6]; - integer_t i, j; - integer_t firstw, lastrw; - float xval, yval; - float* ci; - integer_t isize[2]; - get_dimensions(data, isize); - - assert(state == NULL); - INTERPOLATION_ASSERTS; - - nx = (integer_t)x; - ny = (integer_t)y; - - ci = &coeff[0][0]; - for (j = ny - 2; j <= ny + 3; ++j) { - if (j >= 0 && j < isize[1]) { - for (i = nx - 2; i <= nx + 3; ++i, ++ci) { - if (i < 0) { - *ci = 2.0f * get_pixel(data, 0, j) - get_pixel(data, -i, j); - } else if (i >= isize[0]) { - *ci = 2.0f * get_pixel(data, isize[0]-1, j) - get_pixel(data, 2*isize[0]-2-i, j); + float *value, struct driz_error_t *error UNUSED_PARAM) { + integer_t nx, ny; + const integer_t rowleh = 6; + const integer_t nterms = 6; + float coeff[6][6]; + integer_t i, j; + integer_t firstw, lastrw; + float xval, yval; + float *ci; + integer_t isize[2]; + get_dimensions(data, isize); + + assert(state == NULL); + INTERPOLATION_ASSERTS; + + nx = (integer_t)x; + ny = (integer_t)y; + + ci = &coeff[0][0]; + for (j = ny - 2; j <= ny + 3; ++j) { + if (j >= 0 && j < isize[1]) { + for (i = nx - 2; i <= nx + 3; ++i, ++ci) { + if (i < 0) { + *ci = 2.0f * get_pixel(data, 0, j) - get_pixel(data, -i, j); + } else if (i >= isize[0]) { + *ci = 2.0f * get_pixel(data, isize[0] - 1, j) - + get_pixel(data, 2 * isize[0] - 2 - i, j); + } else { + *ci = get_pixel(data, i, j); + } + } + } else if (j == (ny + 3)) { + for (i = nx - 2; i <= nx + 3; ++i, ++ci) { + if (i < 0) { + *ci = 2.0f * get_pixel(data, 0, isize[1] - 4) - + get_pixel(data, -i, isize[1] - 4); + } else if (i >= isize[0]) { + *ci = 2.0f * get_pixel(data, isize[0] - 1, isize[1] - 4) - + get_pixel(data, 2 * isize[0] - 2 - i, isize[1] - 4); + } else { + *ci = get_pixel(data, i, isize[1] - 4); + } + } } else { - *ci = get_pixel(data, i, j); + ci += 6; } - } - } else if (j == (ny + 3)) { - for (i = nx - 2; i <= nx + 3; ++i, ++ci) { - if (i < 0) { - *ci = 2.0f * get_pixel(data, 0, isize[1]-4) - get_pixel(data, -i, isize[1]-4); - } else if (i >= isize[0]) { - *ci = 2.0f * get_pixel(data, isize[0]-1, isize[1]-4) - get_pixel(data, 2 * isize[0]-2-i, isize[1]-4); - } else { - *ci = get_pixel(data, i, isize[1]-4); - } - } - } else { - ci += 6; } - } - firstw = MAX(0, 2 - ny); - assert(firstw >= 0 && firstw < nterms); + firstw = MAX(0, 2 - ny); + assert(firstw >= 0 && firstw < nterms); - if (firstw > 0) { - for (j = 0; j <= firstw; ++j) { - assert(2*firstw-j >= 0 && 2*firstw-j < nterms); + if (firstw > 0) { + for (j = 0; j <= firstw; ++j) { + assert(2 * firstw - j >= 0 && 2 * firstw - j < nterms); - weighted_sum_vectors(nterms, - &coeff[firstw][0], 2.0, - &coeff[2*firstw-j][0], -1.0, - &coeff[j][0]); + weighted_sum_vectors(nterms, &coeff[firstw][0], 2.0, + &coeff[2 * firstw - j][0], -1.0, &coeff[j][0]); + } } - } - lastrw = MIN(nterms - 1, isize[1] - ny + 1); - assert(lastrw < nterms); + lastrw = MIN(nterms - 1, isize[1] - ny + 1); + assert(lastrw < nterms); + + if (lastrw < nterms - 1) { + for (j = lastrw + 1; j <= nterms - 2; ++j) { + assert(2 * lastrw - j >= 0 && 2 * lastrw - j < nterms); - if (lastrw < nterms - 1) { - for (j = lastrw + 1; j <= nterms - 2; ++j) { - assert(2*lastrw-j >= 0 && 2*lastrw-j < nterms); + weighted_sum_vectors(nterms, &coeff[lastrw][0], 2.0, + &coeff[2 * lastrw - j][0], -1.0, &coeff[j][0]); + } + } else if (lastrw == 2) { + weighted_sum_vectors(nterms, &coeff[2][0], 2.0, &coeff[5][0], -1.0, + &coeff[5][0]); + } else { + assert(2 * lastrw - 5 >= 0 && 2 * lastrw - 5 < nterms); - weighted_sum_vectors(nterms, - &coeff[lastrw][0], 2.0, - &coeff[2*lastrw-j][0], -1.0, - &coeff[j][0]); + weighted_sum_vectors(nterms, &coeff[lastrw][0], 2.0, + &coeff[2 * lastrw - 5][0], -1.0, &coeff[5][0]); } - } else if (lastrw == 2) { - weighted_sum_vectors(nterms, - &coeff[2][0], 2.0, - &coeff[5][0], -1.0, - &coeff[5][0]); - } else { - assert(2*lastrw - 5 >= 0 && 2*lastrw-5 < nterms); - - weighted_sum_vectors(nterms, - &coeff[lastrw][0], 2.0, - &coeff[2*lastrw-5][0], -1.0, - &coeff[5][0]); - } - - xval = 3.0f + (x - (float)nx); - yval = 3.0f + (y - (float)ny); - - ii_bipoly5(&coeff[0][0], rowleh, 0, 1, &xval, &yval, value); - - return 0; + + xval = 3.0f + (x - (float)nx); + yval = 3.0f + (y - (float)ny); + + ii_bipoly5(&coeff[0][0], rowleh, 0, 1, &xval, &yval, value); + + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * (Was: iinisc) */ #define INTERPOLATE_SINC_NCONV 15 static inline_macro int -interpolate_sinc_(PyArrayObject* data, - const integer_t firstt, const integer_t npts, - const float* x /*[npts]*/, const float* y /*[npts]*/, - const float mindx, const float mindy, - const float sinscl, +interpolate_sinc_(PyArrayObject *data, const integer_t firstt, + const integer_t npts, const float *x /*[npts]*/, + const float *y /*[npts]*/, const float mindx, + const float mindy, const float sinscl, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - const integer_t nconv = INTERPOLATE_SINC_NCONV; - const integer_t nsinc = (nconv - 1) / 2; - /* TODO: This is to match Fortan, but is probably technically less precise */ - const float halfpi = 1.5707963267948966192f; /* M_PI / 2.0; */ - const float sconst = powf((halfpi / (float)nsinc), 2.0f); - const float a2 = -0.49670f; - const float a4 = 0.03705f; - float taper[INTERPOLATE_SINC_NCONV]; - float ac[INTERPOLATE_SINC_NCONV], ar[INTERPOLATE_SINC_NCONV]; - float sdx, dx, dy, dxn, dyn, dx2; - float ax, ay, px, py; - float sum, sumx, sumy; - float tmp; - integer_t minj, maxj, offj; - integer_t mink, maxk, offk; - integer_t nx, ny; - integer_t i, j, k, m, index; - integer_t indices[3][4]; - integer_t isize[2]; - get_dimensions(data, isize); - - assert(x); - assert(y); - assert(value); - assert(error); - - if ((nsinc % 2) == 0) { - sdx = 1.0; - for (j = -nsinc; j <= nsinc; ++j) { - assert(j + nsinc >= 0 && j + nsinc < INTERPOLATE_SINC_NCONV); - - taper[j + nsinc] = 1.0; + float *value, struct driz_error_t *error UNUSED_PARAM) { + const integer_t nconv = INTERPOLATE_SINC_NCONV; + const integer_t nsinc = (nconv - 1) / 2; + /* TODO: This is to match Fortan, but is probably technically less precise + */ + const float halfpi = 1.5707963267948966192f; /* M_PI / 2.0; */ + const float sconst = powf((halfpi / (float)nsinc), 2.0f); + const float a2 = -0.49670f; + const float a4 = 0.03705f; + float taper[INTERPOLATE_SINC_NCONV]; + float ac[INTERPOLATE_SINC_NCONV], ar[INTERPOLATE_SINC_NCONV]; + float sdx, dx, dy, dxn, dyn, dx2; + float ax, ay, px, py; + float sum, sumx, sumy; + float tmp; + integer_t minj, maxj, offj; + integer_t mink, maxk, offk; + integer_t nx, ny; + integer_t i, j, k, m, index; + integer_t indices[3][4]; + integer_t isize[2]; + get_dimensions(data, isize); + + assert(x); + assert(y); + assert(value); + assert(error); + + if ((nsinc % 2) == 0) { + sdx = 1.0; + for (j = -nsinc; j <= nsinc; ++j) { + assert(j + nsinc >= 0 && j + nsinc < INTERPOLATE_SINC_NCONV); + + taper[j + nsinc] = 1.0; + } + } else { + sdx = -1.0; + errno = 0; + for (j = -nsinc; j <= nsinc; ++j) { + assert(j + nsinc >= 0 && j + nsinc < INTERPOLATE_SINC_NCONV); + + dx2 = sconst * (float)j * (float)j; + tmp = powf(1.0f + a2 * dx2 + a4 * dx2 * dx2, 2.0); + if (errno != 0) { + driz_error_set_message(error, "pow failed"); + return 1; + } + taper[j + nsinc] = sdx * tmp; + + sdx = -sdx; + } } - } else { - sdx = -1.0; - errno = 0; - for (j = -nsinc; j <= nsinc; ++j) { - assert(j + nsinc >= 0 && j + nsinc < INTERPOLATE_SINC_NCONV); - - dx2 = sconst * (float)j * (float)j; - tmp = powf(1.0f + a2*dx2 + a4*dx2*dx2, 2.0); - if (errno != 0) { - driz_error_set_message(error, "pow failed"); - return 1; - } - taper[j + nsinc] = sdx * tmp; - sdx = -sdx; - } - } - - for (i = 0; i < npts; ++i) { - nx = fortran_round(x[i]); - ny = fortran_round(y[i]); - if (nx < 0 || nx >= isize[0] || ny < 0 || ny >= isize[1]) { - value[i] = 0.0; - continue; - } + for (i = 0; i < npts; ++i) { + nx = fortran_round(x[i]); + ny = fortran_round(y[i]); + if (nx < 0 || nx >= isize[0] || ny < 0 || ny >= isize[1]) { + value[i] = 0.0; + continue; + } - dx = (x[i] - (float)nx) * sinscl; - dy = (y[i] - (float)ny) * sinscl; + dx = (x[i] - (float)nx) * sinscl; + dy = (y[i] - (float)ny) * sinscl; - if (fabsf(dx) < mindx && fabsf(dy) < mindy) { - index = firstt + (ny - 1) * isize[0] + nx - 1; /* TODO: Base check */ - value[i] = get_pixel_at_pos(data, index); - continue; - } + if (fabsf(dx) < mindx && fabsf(dy) < mindy) { + index = + firstt + (ny - 1) * isize[0] + nx - 1; /* TODO: Base check */ + value[i] = get_pixel_at_pos(data, index); + continue; + } - dxn = 1.0f + (float)nsinc + dx; - dyn = 1.0f + (float)nsinc + dy; - sumx = 0.0f; - sumy = 0.0f; - for (j = 0; j < nconv; ++j) { - /* TODO: These out of range indices also seem to be in Fortran... */ - ax = dxn - (float)j - 1; /* TODO: Base check */ - ay = dyn - (float)j - 1; /* TODO: Base check */ - assert(ax != 0.0); - assert(ay != 0.0); - - if (ax == 0.0) { - px = 1.0; - } else if (dx == 0.0) { - px = 0.0; - } else { - px = taper[j - 1] / ax; - } - - if (ay == 0.0) { - py = 1.0; - } else if (dy == 0.0) { - py = 0.0; - } else { - py = taper[j - 1] / ay; - } - - /* TODO: These out of range indices also seem to be in Fortran... */ - ac[j - 1] = px; - ar[j - 1] = py; - sumx += px; - sumy += py; - } + dxn = 1.0f + (float)nsinc + dx; + dyn = 1.0f + (float)nsinc + dy; + sumx = 0.0f; + sumy = 0.0f; + for (j = 0; j < nconv; ++j) { + /* TODO: These out of range indices also seem to be in Fortran... */ + ax = dxn - (float)j - 1; /* TODO: Base check */ + ay = dyn - (float)j - 1; /* TODO: Base check */ + assert(ax != 0.0); + assert(ay != 0.0); + + if (ax == 0.0) { + px = 1.0; + } else if (dx == 0.0) { + px = 0.0; + } else { + px = taper[j - 1] / ax; + } + + if (ay == 0.0) { + py = 1.0; + } else if (dy == 0.0) { + py = 0.0; + } else { + py = taper[j - 1] / ay; + } + + /* TODO: These out of range indices also seem to be in Fortran... */ + ac[j - 1] = px; + ar[j - 1] = py; + sumx += px; + sumy += py; + } - /* Compute the limits of the convolution */ - minj = MAX(0, ny - nsinc - 1); /* TODO: Bases check */ - maxj = MIN(isize[1], ny + nsinc); /* TODO: Bases check */ - offj = nsinc - ny; /* TODO: Bases check */ + /* Compute the limits of the convolution */ + minj = MAX(0, ny - nsinc - 1); /* TODO: Bases check */ + maxj = MIN(isize[1], ny + nsinc); /* TODO: Bases check */ + offj = nsinc - ny; /* TODO: Bases check */ - mink = MAX(0, nx - nsinc - 1); /* TODO: Bases check */ - maxk = MIN(isize[0], nx + nsinc); /* TODO: Bases check */ - offk = nsinc - nx; /* TODO: Bases check */ + mink = MAX(0, nx - nsinc - 1); /* TODO: Bases check */ + maxk = MIN(isize[0], nx + nsinc); /* TODO: Bases check */ + offk = nsinc - nx; /* TODO: Bases check */ - value[i] = 0.0; + value[i] = 0.0; - indices[0][0] = ny - nsinc; - indices[0][1] = minj - 1; - indices[0][2] = firstt; - indices[0][3] = 0; + indices[0][0] = ny - nsinc; + indices[0][1] = minj - 1; + indices[0][2] = firstt; + indices[0][3] = 0; - indices[1][0] = minj; - indices[1][1] = maxj; - indices[1][2] = firstt; - indices[1][3] = isize[0]; + indices[1][0] = minj; + indices[1][1] = maxj; + indices[1][2] = firstt; + indices[1][3] = isize[0]; - indices[2][0] = maxj + 1; - indices[2][1] = ny + nsinc; - indices[2][2] = firstt + (isize[1] - 1) * isize[0]; - indices[2][3] = 0; + indices[2][0] = maxj + 1; + indices[2][1] = ny + nsinc; + indices[2][2] = firstt + (isize[1] - 1) * isize[0]; + indices[2][3] = 0; - for (m = 0; m < 3; ++m) { - for (j = indices[m][0]; j <= indices[m][1]; ++j) { - sum = 0.0; - index = indices[m][2] + j * indices[m][3]; - assert(index >= 0 && index < isize[0]*isize[1] - 1); - assert(index+isize[0] >= 0 && index+isize[0] < isize[0]*isize[1]); + for (m = 0; m < 3; ++m) { + for (j = indices[m][0]; j <= indices[m][1]; ++j) { + sum = 0.0; + index = indices[m][2] + j * indices[m][3]; + assert(index >= 0 && index < isize[0] * isize[1] - 1); + assert(index + isize[0] >= 0 && + index + isize[0] < isize[0] * isize[1]); - for (k = nx - nsinc; k < mink - 1; ++k) { /* TODO: Bases check */ - assert(k+offk >= 0 && k+offk < INTERPOLATE_SINC_NCONV); + for (k = nx - nsinc; k < mink - 1; + ++k) { /* TODO: Bases check */ + assert(k + offk >= 0 && k + offk < INTERPOLATE_SINC_NCONV); - sum += ac[k+offk] * get_pixel_at_pos(data, index+1); - } + sum += ac[k + offk] * get_pixel_at_pos(data, index + 1); + } - for (k = mink; k <= maxk; ++k) { /* TODO: Bases check */ - assert(k+offk >= 0 && k+offk < INTERPOLATE_SINC_NCONV); - assert(index+k >= 0 && index+k < isize[0]*isize[1]); + for (k = mink; k <= maxk; ++k) { /* TODO: Bases check */ + assert(k + offk >= 0 && k + offk < INTERPOLATE_SINC_NCONV); + assert(index + k >= 0 && index + k < isize[0] * isize[1]); - sum += ac[k+offk] * get_pixel_at_pos(data, index+k); - } + sum += ac[k + offk] * get_pixel_at_pos(data, index + k); + } - for (k = maxk + 1; k <= nx + nsinc; ++k) { - assert(k+offk >= 0 && k+offk < INTERPOLATE_SINC_NCONV); + for (k = maxk + 1; k <= nx + nsinc; ++k) { + assert(k + offk >= 0 && k + offk < INTERPOLATE_SINC_NCONV); - sum += ac[k+offk] * get_pixel_at_pos(data, index+isize[0]); - } + sum += + ac[k + offk] * get_pixel_at_pos(data, index + isize[0]); + } - assert(j + offj >= 0 && j + offj < INTERPOLATE_SINC_NCONV); + assert(j + offj >= 0 && j + offj < INTERPOLATE_SINC_NCONV); - value[i] += ar[j + offj] * sum; - } - } + value[i] += ar[j + offj] * sum; + } + } - assert(sumx != 0.0); - assert(sumy != 0.0); + assert(sumx != 0.0); + assert(sumy != 0.0); - value[i] = value[i] / sumx / sumy; - } + value[i] = value[i] / sumx / sumy; + } - return 0; + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform sinc interpolation. * * state: A pointer to a \a sinc_param_t object @@ -767,24 +741,22 @@ interpolate_sinc_(PyArrayObject* data, */ static int -interpolate_sinc(const void* state, - PyArrayObject* data, - const float x, const float y, +interpolate_sinc(const void *state, PyArrayObject *data, const float x, + const float y, /* Output parameters */ - float* value, - struct driz_error_t* error) { - const struct sinc_param_t* param = (const struct sinc_param_t*)state; - integer_t isize[2]; - get_dimensions(data, isize); + float *value, struct driz_error_t *error) { + const struct sinc_param_t *param = (const struct sinc_param_t *)state; + integer_t isize[2]; + get_dimensions(data, isize); - assert(state); - INTERPOLATION_ASSERTS; + assert(state); + INTERPOLATION_ASSERTS; - return interpolate_sinc_(data, 0, 1, &x, &y, 0.001f, 0.001f, - param->sinscl, value, error); + return interpolate_sinc_(data, 0, 1, &x, &y, 0.001f, 0.001f, param->sinscl, + value, error); } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Perform Lanczos interpolation. * * state: A pointer to a \a lanczos_param_t object @@ -796,203 +768,204 @@ interpolate_sinc(const void* state, */ static int -interpolate_lanczos(const void* state, - PyArrayObject* data, - const float x, const float y, +interpolate_lanczos(const void *state, PyArrayObject *data, const float x, + const float y, /* Output parameters */ - float* value, - struct driz_error_t* error UNUSED_PARAM) { - integer_t ixs, iys, ixe, iye; - integer_t xoff, yoff; - float luty, sum; - integer_t nbox; - integer_t i, j; - const struct lanczos_param_t* lanczos = (const struct lanczos_param_t*)state; - integer_t isize[2]; - get_dimensions(data, isize); - - assert(state); - INTERPOLATION_ASSERTS; - - nbox = lanczos->nbox; - - /* First check for being close to the edge and, if so, return the - missing value */ - ixs = (integer_t)(x) - nbox; - ixe = (integer_t)(x) + nbox; - iys = (integer_t)(y) - nbox; - iye = (integer_t)(y) + nbox; - if (ixs < 0 || ixe >= isize[0] || - iys < 0 || iye >= isize[1]) { - *value = lanczos->misval; - return 0; - } + float *value, struct driz_error_t *error UNUSED_PARAM) { + integer_t ixs, iys, ixe, iye; + integer_t xoff, yoff; + float luty, sum; + integer_t nbox; + integer_t i, j; + const struct lanczos_param_t *lanczos = + (const struct lanczos_param_t *)state; + integer_t isize[2]; + get_dimensions(data, isize); + + assert(state); + INTERPOLATION_ASSERTS; + + nbox = lanczos->nbox; + + /* First check for being close to the edge and, if so, return the + missing value */ + ixs = (integer_t)(x)-nbox; + ixe = (integer_t)(x) + nbox; + iys = (integer_t)(y)-nbox; + iye = (integer_t)(y) + nbox; + if (ixs < 0 || ixe >= isize[0] || iys < 0 || iye >= isize[1]) { + *value = lanczos->misval; + return 0; + } - /* Don't divide-by-zero errors */ - assert(lanczos->space != 0.0); + /* Don't divide-by-zero errors */ + assert(lanczos->space != 0.0); - /* Loop over the box, which is assumed to be scaled appropriately */ - sum = 0.0; - for (j = iys; j <= iye; ++j) { - yoff = (integer_t)(fabs((y - (float)j) / lanczos->space)); - assert(yoff >= 0 && yoff < lanczos->nlut); + /* Loop over the box, which is assumed to be scaled appropriately */ + sum = 0.0; + for (j = iys; j <= iye; ++j) { + yoff = (integer_t)(fabs((y - (float)j) / lanczos->space)); + assert(yoff >= 0 && yoff < lanczos->nlut); - luty = lanczos->lut[yoff]; - for (i = ixs; i <= ixe; ++i) { - xoff = (integer_t)(fabs((x - (float)i) / lanczos->space)); - assert(xoff >= 0 && xoff < lanczos->nlut); + luty = lanczos->lut[yoff]; + for (i = ixs; i <= ixe; ++i) { + xoff = (integer_t)(fabs((x - (float)i) / lanczos->space)); + assert(xoff >= 0 && xoff < lanczos->nlut); - sum += get_pixel(data, i, j) * lanczos->lut[xoff] * luty; + sum += get_pixel(data, i, j) * lanczos->lut[xoff] * luty; + } } - } - *value = sum; - return 0; + *value = sum; + return 0; } -/** -------------------------------------------------------------------------------------------------- - * A mapping from e_interp_t enumeration values to function pointers that actually - * perform the interpolation. NULL elements will raise an "unimplemented" error. +/** --------------------------------------------------------------------------- + * A mapping from e_interp_t enumeration values to function pointers that + * actually perform the interpolation. NULL elements will raise an + * "unimplemented" error. */ -interp_function* interp_function_map[interp_LAST] = { - &interpolate_nearest_neighbor, - &interpolate_bilinear, - &interpolate_poly3, - &interpolate_poly5, - NULL, - &interpolate_sinc, - &interpolate_sinc, - &interpolate_lanczos, - &interpolate_lanczos -}; - -/** -------------------------------------------------------------------------------------------------- +interp_function *interp_function_map[interp_LAST] = { + &interpolate_nearest_neighbor, + &interpolate_bilinear, + &interpolate_poly3, + &interpolate_poly5, + NULL, + &interpolate_sinc, + &interpolate_sinc, + &interpolate_lanczos, + &interpolate_lanczos}; + +/** --------------------------------------------------------------------------- * Interpolate grid of pixels onto new grid of different size. * * p: structure containing options, input, and output */ int -doblot(struct driz_param_t* p) { - - const size_t nlut = 2048; - const float space = 0.01; - integer_t isize[2], osize[2]; - float scale2, xo, yo, v; - integer_t i, j; - interp_function* interpolate; - struct sinc_param_t sinc; - struct lanczos_param_t lanczos; - void* state = NULL; - - driz_log_message("starting doblot"); - get_dimensions(p->data, isize); - get_dimensions(p->output_data, osize); - - /* Select interpolation function */ - assert(p->interpolation >= 0 && p->interpolation < interp_LAST); - interpolate = interp_function_map[p->interpolation]; - if (interpolate == NULL) { - driz_error_set_message(p->error, "Requested interpolation type not implemented."); - goto doblot_exit_; - } - - lanczos.lut = NULL; - - /* Some interpolation functions need some pre-calculated state */ - if (p->interpolation == interp_lanczos3 || p->interpolation == interp_lanczos5) { - - if ((lanczos.lut = (float*)malloc(nlut * sizeof(float))) == NULL) { - driz_error_set_message(p->error, "Out of memory"); - goto doblot_exit_; +doblot(struct driz_param_t *p) { + const size_t nlut = 2048; + const float space = 0.01; + integer_t isize[2], osize[2]; + float scale2, xo, yo, v; + integer_t i, j; + interp_function *interpolate; + struct sinc_param_t sinc; + struct lanczos_param_t lanczos; + void *state = NULL; + + driz_log_message("starting doblot"); + get_dimensions(p->data, isize); + get_dimensions(p->output_data, osize); + + /* Select interpolation function */ + assert(p->interpolation >= 0 && p->interpolation < interp_LAST); + interpolate = interp_function_map[p->interpolation]; + if (interpolate == NULL) { + driz_error_set_message(p->error, + "Requested interpolation type not implemented."); + goto doblot_exit_; } - create_lanczos_lut(p->interpolation == interp_lanczos3 ? 3 : 5, - nlut, space, lanczos.lut); - - lanczos.nbox = (integer_t)(3.0 / p->kscale); - lanczos.nlut = nlut; - lanczos.space = space; - lanczos.misval = p->misval; - - state = &lanczos; - - } else if (p->interpolation == interp_sinc || p->interpolation == interp_lsinc) { - sinc.sinscl = p->sinscl; - state = &sinc; - - } /* Otherwise state is NULL */ - - /* In the WCS case, we can't use the scale to calculate the Jacobian, - so we need to do it. - - Note that we use the center of the image, rather than the reference pixel - as the reference here. + lanczos.lut = NULL; - This is taken from dobox, except for the inversion of the image order. - - This section applies in WBLOT mode and now contains the addition - correction to separate the distortion-induced scale change. - */ - - /* Recalculate the area scaling factor */ - scale2 = p->scale * p->scale; - v = 1.0; - - for (j = 0; j < osize[1]; ++j) { - - /* Loop through the output positions and do the interpolation */ - for (i = 0; i < osize[0]; ++i) { - if (oob_pixel(p->pixmap, i, j)) { - driz_error_format_message(p->error, "OOB in pixmap[%d,%d]", i, j); - return 1; - } else { - xo = get_pixmap(p->pixmap, i, j)[0]; - yo = get_pixmap(p->pixmap, i, j)[1]; - } - - if (npy_isnan(xo) || npy_isnan(yo)) { - driz_error_format_message(p->error, "NaN in pixmap[%d,%d]", i, j); - return 1; - } - - /* Check it is on the input image */ - if (xo >= 0.0 && xo < (float)isize[0] && - yo >= 0.0 && yo < (float)isize[1]) { - - double value; - - /* Check for look-up-table interpolation */ - if (interpolate(state, p->data, xo, yo, &v, p->error)) { - goto doblot_exit_; - } - - value = v * p->ef / scale2; - if (oob_pixel(p->output_data, i, j)) { - driz_error_format_message(p->error, "OOB in output_data[%d,%d]", i, j); - return 1; - } else { - set_pixel(p->output_data, i, j, value); + /* Some interpolation functions need some pre-calculated state */ + if (p->interpolation == interp_lanczos3 || + p->interpolation == interp_lanczos5) { + if ((lanczos.lut = (float *)malloc(nlut * sizeof(float))) == NULL) { + driz_error_set_message(p->error, "Out of memory"); + goto doblot_exit_; } - } else { - /* If there is nothing for us then set the output to missing C - value flag */ - if (oob_pixel(p->output_data, i, j)) { - driz_error_format_message(p->error, "OOB in output_data[%d,%d]", i, j); - return 1; - } else { - set_pixel(p->output_data, i, j, p->misval); - p->nmiss++; + create_lanczos_lut(p->interpolation == interp_lanczos3 ? 3 : 5, nlut, + space, lanczos.lut); + + lanczos.nbox = (integer_t)(3.0 / p->kscale); + lanczos.nlut = nlut; + lanczos.space = space; + lanczos.misval = p->misval; + + state = &lanczos; + + } else if (p->interpolation == interp_sinc || + p->interpolation == interp_lsinc) { + sinc.sinscl = p->sinscl; + state = &sinc; + + } /* Otherwise state is NULL */ + + /* In the WCS case, we can't use the scale to calculate the Jacobian, + so we need to do it. + + Note that we use the center of the image, rather than the reference pixel + as the reference here. + + This is taken from dobox, except for the inversion of the image order. + + This section applies in WBLOT mode and now contains the addition + correction to separate the distortion-induced scale change. + */ + + /* Recalculate the area scaling factor */ + scale2 = p->scale * p->scale; + v = 1.0; + + for (j = 0; j < osize[1]; ++j) { + /* Loop through the output positions and do the interpolation */ + for (i = 0; i < osize[0]; ++i) { + if (oob_pixel(p->pixmap, i, j)) { + driz_error_format_message(p->error, "OOB in pixmap[%d,%d]", i, + j); + return 1; + } else { + xo = get_pixmap(p->pixmap, i, j)[0]; + yo = get_pixmap(p->pixmap, i, j)[1]; + } + + if (npy_isnan(xo) || npy_isnan(yo)) { + driz_error_format_message(p->error, "NaN in pixmap[%d,%d]", i, + j); + return 1; + } + + /* Check it is on the input image */ + if (xo >= 0.0 && xo < (float)isize[0] && yo >= 0.0 && + yo < (float)isize[1]) { + double value; + + /* Check for look-up-table interpolation */ + if (interpolate(state, p->data, xo, yo, &v, p->error)) { + goto doblot_exit_; + } + + value = v * p->ef / scale2; + if (oob_pixel(p->output_data, i, j)) { + driz_error_format_message( + p->error, "OOB in output_data[%d,%d]", i, j); + return 1; + } else { + set_pixel(p->output_data, i, j, value); + } + + } else { + /* If there is nothing for us then set the output to missing C + value flag */ + if (oob_pixel(p->output_data, i, j)) { + driz_error_format_message( + p->error, "OOB in output_data[%d,%d]", i, j); + return 1; + } else { + set_pixel(p->output_data, i, j, p->misval); + p->nmiss++; + } + } } - } } - } - doblot_exit_: - driz_log_message("ending doblot"); - if (lanczos.lut) free(lanczos.lut); +doblot_exit_: + driz_log_message("ending doblot"); + if (lanczos.lut) free(lanczos.lut); - return driz_error_is_set(p->error); + return driz_error_is_set(p->error); } diff --git a/src/cdrizzleblot.h b/src/cdrizzleblot.h index d07c4b41..58b9fb4d 100644 --- a/src/cdrizzleblot.h +++ b/src/cdrizzleblot.h @@ -10,7 +10,6 @@ * @return Non-zero if an error occurred. */ -int -doblot(struct driz_param_t* p); +int doblot(struct driz_param_t *p); #endif /* CDRIZZLEBLOT_H */ diff --git a/src/cdrizzlebox.c b/src/cdrizzlebox.c index f6b401d1..b29b8ef6 100644 --- a/src/cdrizzlebox.c +++ b/src/cdrizzlebox.c @@ -7,12 +7,11 @@ #include "cdrizzleutil.h" #include -#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ +#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ #include #include - -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * Update the flux and counts in the output image using a weighted average * * p: structure containing options, input, and output @@ -24,45 +23,48 @@ */ inline_macro static int -update_data(struct driz_param_t* p, const integer_t ii, const integer_t jj, +update_data(struct driz_param_t *p, const integer_t ii, const integer_t jj, const float d, const float vc, const float dow) { + double vc_plus_dow; - double vc_plus_dow; + if (dow == 0.0f) return 0; - if (dow == 0.0f) return 0; + vc_plus_dow = vc + dow; - vc_plus_dow = vc + dow; + if (vc == 0.0f) { + if (oob_pixel(p->output_data, ii, jj)) { + driz_error_format_message(p->error, "OOB in output_data[%d,%d]", ii, + jj); + return 1; + } else { + set_pixel(p->output_data, ii, jj, d); + } - if (vc == 0.0f) { - if (oob_pixel(p->output_data, ii, jj)) { - driz_error_format_message(p->error, "OOB in output_data[%d,%d]", ii, jj); - return 1; } else { - set_pixel(p->output_data, ii, jj, d); + if (oob_pixel(p->output_data, ii, jj)) { + driz_error_format_message(p->error, "OOB in output_data[%d,%d]", ii, + jj); + return 1; + } else { + double value; + value = (get_pixel(p->output_data, ii, jj) * vc + dow * d) / + (vc_plus_dow); + set_pixel(p->output_data, ii, jj, value); + } } - } else { - if (oob_pixel(p->output_data, ii, jj)) { - driz_error_format_message(p->error, "OOB in output_data[%d,%d]", ii, jj); - return 1; + if (oob_pixel(p->output_counts, ii, jj)) { + driz_error_format_message(p->error, "OOB in output_counts[%d,%d]", ii, + jj); + return 1; } else { - double value; - value = (get_pixel(p->output_data, ii, jj) * vc + dow * d) / (vc_plus_dow); - set_pixel(p->output_data, ii, jj, value); + set_pixel(p->output_counts, ii, jj, vc_plus_dow); } - } - - if (oob_pixel(p->output_counts, ii, jj)) { - driz_error_format_message(p->error, "OOB in output_counts[%d,%d]", ii, jj); - return 1; - } else { - set_pixel(p->output_counts, ii, jj, vc_plus_dow); - } - return 0; + return 0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * The bit value, trimmed to the appropriate range * * uuid: the id of the input image @@ -80,172 +82,180 @@ compute_bit_value(integer_t uuid) { return bv; } -/** -------------------------------------------------------------------------------------------------- - * Compute area of box overlap. Calculate the area common to input clockwise polygon x(n), y(n) with - * square (is, js) to (is+1, js+1). This version is for a quadrilateral. Used by do_square_kernel. +/** --------------------------------------------------------------------------- + * Compute area of box overlap. Calculate the area common to input clockwise + * polygon x(n), y(n) with square (is, js) to (is+1, js+1). This version is for + * a quadrilateral. Used by do_square_kernel. * * is: x coordinate of a pixel on the output image * js: y coordinate of a pixel on the output image - * x: x coordinates of endpoints of quadrilateral containing flux of input pixel - * y: y coordinates of endpoints of quadrilateral containing flux of input pixel + * x: x coordinates of endpoints of quadrilateral containing flux of input + * pixel y: y coordinates of endpoints of quadrilateral containing flux of + * input pixel */ double compute_area(double is, double js, const double x[4], const double y[4]) { - int ipoint, jpoint, idim, jdim, iside, outside, count; - int positive[2]; - double area, width; - double midpoint[2], delta[2]; - double border[2][2], segment[2][2]; - - /* The area for a qadrilateral clipped to a square of unit length whose sides are - * aligned with the axes. The area is computed by computing the area under each - * line segment clipped to the boundary of three sides of the sqaure. Since the - * computed width is positive for two of the sides and negative for the other two, - * we subtract the area outside the quadrilateral without any extra code. - */ - area = 0.0; - - border[0][0] = is - 0.5; - border[0][1] = js - 0.5; - border[1][0] = is + 0.5; - border[1][1] = js + 0.5; - - for (ipoint = 0; ipoint < 4; ++ ipoint) { - jpoint = (ipoint + 1) & 03; /* Next point in cyclical order */ - - segment[0][0] = x[ipoint]; - segment[0][1] = y[ipoint]; - segment[1][0] = x[jpoint]; - segment[1][1] = y[jpoint]; - - /* Compute the endpoints of the line segment that - * lie inside the border (possibly the whole segment) + int ipoint, jpoint, idim, jdim, iside, outside, count; + int positive[2]; + double area, width; + double midpoint[2], delta[2]; + double border[2][2], segment[2][2]; + + /* The area for a qadrilateral clipped to a square of unit length whose + * sides are aligned with the axes. The area is computed by computing the + * area under each line segment clipped to the boundary of three sides of + * the sqaure. Since the computed width is positive for two of the sides and + * negative for the other two, we subtract the area outside the + * quadrilateral without any extra code. */ + area = 0.0; - for (idim = 0, count = 3; idim < 2; ++ idim) { - for (iside = 0; iside < 2; ++ iside, -- count) { + border[0][0] = is - 0.5; + border[0][1] = js - 0.5; + border[1][0] = is + 0.5; + border[1][1] = js + 0.5; - delta[0] = segment[0][idim] - border[iside][idim]; - delta[1] = segment[1][idim] - border[iside][idim]; + for (ipoint = 0; ipoint < 4; ++ipoint) { + jpoint = (ipoint + 1) & 03; /* Next point in cyclical order */ - positive[0] = delta[0] > 0.0; - positive[1] = delta[1] > 0.0; + segment[0][0] = x[ipoint]; + segment[0][1] = y[ipoint]; + segment[1][0] = x[jpoint]; + segment[1][1] = y[jpoint]; - /* If both deltas have the same signe there is no baundary crossing + /* Compute the endpoints of the line segment that + * lie inside the border (possibly the whole segment) */ - if (positive[0] == positive[1]) { - /* A diagram will convince that you decide a point is - * inside or outside the boundary by the following test - */ - if (positive[0] == iside) { - /* Segment is entirely outside the boundary */ - if (count == 0) { - /* Implicitly multiplied by 1.0, the square height */ - width = segment[1][0] - segment[0][0]; - area += width; - } else { - goto _nextsegment; - } - } else { - /* Segment entirely within the boundary */ - if (count == 0) { - /* Use the trapezoid formula to compute the area under the - * segment. Delta is the distance to the top of the square - * and is negative or zero for the segment inside the square - */ - width = segment[1][0] - segment[0][0]; - area += 0.5 * width * ((1.0 + delta[0]) + (1.0 + delta[1])); - } - } + for (idim = 0, count = 3; idim < 2; ++idim) { + for (iside = 0; iside < 2; ++iside, --count) { + delta[0] = segment[0][idim] - border[iside][idim]; + delta[1] = segment[1][idim] - border[iside][idim]; + + positive[0] = delta[0] > 0.0; + positive[1] = delta[1] > 0.0; + + /* If both deltas have the same signe there is no baundary + * crossing + */ + if (positive[0] == positive[1]) { + /* A diagram will convince that you decide a point is + * inside or outside the boundary by the following test + */ + if (positive[0] == iside) { + /* Segment is entirely outside the boundary */ + if (count == 0) { + /* Implicitly multiplied by 1.0, the square height + */ + width = segment[1][0] - segment[0][0]; + area += width; + } else { + goto _nextsegment; + } - } else { - /* If ends of the line segment are on opposite sides of the - * boundary, calculate midpoint, the point of intersection - */ - outside = positive[iside]; - jdim = (idim + 1) & 01; /* the other dimension */ - - midpoint[idim] = border[iside][idim]; - - midpoint[jdim] = - (delta[1] * segment[0][jdim] - delta[0] * segment[1][jdim]) / - (delta[1] - delta[0]); - - if (count == 0) { - /* If a segment cross the boundary the formula for its area - * is a combination of the formulas for segments entirely - * inside and outside. - */ - if (outside == 0) { - width = midpoint[0] - segment[0][0]; - area += width; - width = segment[1][0] - midpoint[0]; - /* Delta[0] is at the crossing point and thus zero */ - area += 0.5 * width * (1.0 + (1.0 + delta[1])); - } else { - width = segment[1][0] - midpoint[0]; - area += width; - width = midpoint[0] - segment[0][0]; - /* Delta[1] is at the crossing point and thus zero */ - area += 0.5 * width * ((1.0 + delta[0]) + 1.0); - } + } else { + /* Segment entirely within the boundary */ + if (count == 0) { + /* Use the trapezoid formula to compute the area + * under the segment. Delta is the distance to the + * top of the square and is negative or zero for the + * segment inside the square + */ + width = segment[1][0] - segment[0][0]; + area += 0.5 * width * + ((1.0 + delta[0]) + (1.0 + delta[1])); + } + } + + } else { + /* If ends of the line segment are on opposite sides of the + * boundary, calculate midpoint, the point of intersection + */ + outside = positive[iside]; + jdim = (idim + 1) & 01; /* the other dimension */ + + midpoint[idim] = border[iside][idim]; + + midpoint[jdim] = (delta[1] * segment[0][jdim] - + delta[0] * segment[1][jdim]) / + (delta[1] - delta[0]); + + if (count == 0) { + /* If a segment cross the boundary the formula for its + * area is a combination of the formulas for segments + * entirely inside and outside. + */ + if (outside == 0) { + width = midpoint[0] - segment[0][0]; + area += width; + width = segment[1][0] - midpoint[0]; + /* Delta[0] is at the crossing point and thus zero + */ + area += 0.5 * width * (1.0 + (1.0 + delta[1])); + } else { + width = segment[1][0] - midpoint[0]; + area += width; + width = midpoint[0] - segment[0][0]; + /* Delta[1] is at the crossing point and thus zero + */ + area += 0.5 * width * ((1.0 + delta[0]) + 1.0); + } - } else { - /* Clip segment against each boundary except the last */ - segment[outside][0] = midpoint[0]; - segment[outside][1] = midpoint[1]; - } + } else { + /* Clip segment against each boundary except the last */ + segment[outside][0] = midpoint[0]; + segment[outside][1] = midpoint[1]; + } + } + } } - } - } - _nextsegment: continue; - } + _nextsegment: + continue; + } - return fabs(area); + return fabs(area); } -/** -------------------------------------------------------------------------------------------------- - * Calculate overlap between an arbitrary rectangle, aligned with the axes, and a pixel. - * This is a simplified version of the compute_area, only valid if axes are nearly aligned. - * Used by do_kernel_turbo. +/** --------------------------------------------------------------------------- + * Calculate overlap between an arbitrary rectangle, aligned with the axes, and + * a pixel. This is a simplified version of the compute_area, only valid if axes + * are nearly aligned. Used by do_kernel_turbo. * * i: the x coordinate of a pixel on the output image * j: the y coordinate of a pixel on the output image - * xmin: the x coordinate of the lower edge of rectangle containing flux of input pixel - * xmax: the x coordinate of the upper edge of rectangle containing flux of input pixel - * ymin: the y coordinate of the lower edge of rectangle containing flux of input pixel - * ymax: the y coordinate of the upper edge of rectangle containing flux of input pixel + * xmin: the x coordinate of the lower edge of rectangle containing flux of + * input pixel xmax: the x coordinate of the upper edge of rectangle containing + * flux of input pixel ymin: the y coordinate of the lower edge of rectangle + * containing flux of input pixel ymax: the y coordinate of the upper edge of + * rectangle containing flux of input pixel */ static inline_macro double -over(const integer_t i, const integer_t j, - const double xmin, const double xmax, +over(const integer_t i, const integer_t j, const double xmin, const double xmax, const double ymin, const double ymax) { double dx, dy; assert(xmin <= xmax); assert(ymin <= ymax); - dx = MIN(xmax, (double)(i) + 0.5) - MAX(xmin, (double)(i) - 0.5); - dy = MIN(ymax, (double)(j) + 0.5) - MAX(ymin, (double)(j) - 0.5); + dx = MIN(xmax, (double)(i) + 0.5) - MAX(xmin, (double)(i)-0.5); + dy = MIN(ymax, (double)(j) + 0.5) - MAX(ymin, (double)(j)-0.5); - if (dx > 0.0 && dy > 0.0) - return dx*dy; + if (dx > 0.0 && dy > 0.0) return dx * dy; return 0.0; } -/** -------------------------------------------------------------------------------------------------- +/** --------------------------------------------------------------------------- * The kernel assumes all the flux in an input pixel is at the center * * p: structure containing options, input, and output */ static int -do_kernel_point(struct driz_param_t* p) { +do_kernel_point(struct driz_param_t *p) { struct scanner s; integer_t i, j, ii, jj; integer_t osize[2]; @@ -272,8 +282,8 @@ do_kernel_point(struct driz_param_t* p) { p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); break; } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) + // pixel centered on y is outside of scanner's limits or image [0, + // height - 1] OR: limits (x1, x2) are equal (line width is 0) p->nmiss += (p->xmax - p->xmin); ++p->nskip; continue; @@ -286,7 +296,7 @@ do_kernel_point(struct driz_param_t* p) { double ox, oy; if (map_pixel(p->pixmap, i, j, &ox, &oy)) { - ++ p->nmiss; + ++p->nmiss; } else { ii = fortran_round(ox); @@ -294,7 +304,7 @@ do_kernel_point(struct driz_param_t* p) { /* Check it is on the output image */ if (ii < 0 || ii >= osize[0] || jj < 0 || jj >= osize[1]) { - ++ p->nmiss; + ++p->nmiss; } else { vc = get_pixel(p->output_counts, ii, jj); @@ -302,8 +312,9 @@ do_kernel_point(struct driz_param_t* p) { /* Allow for stretching because of scale change */ d = get_pixel(p->data, i, j) * scale2; - /* Scale the weighting mask by the scale factor. Note that we - DON'T scale by the Jacobian as it hasn't been calculated */ + /* Scale the weighting mask by the scale factor. Note that + we DON'T scale by the Jacobian as it hasn't been + calculated */ if (p->weights) { dow = get_pixel(p->weights, i, j) * p->weight_scale; } else { @@ -327,21 +338,21 @@ do_kernel_point(struct driz_param_t* p) { return 0; } - -/** -------------------------------------------------------------------------------------------------- - * This kernel assumes the flux is distributed acrass a gaussian around the center of an input pixel +/** --------------------------------------------------------------------------- + * This kernel assumes the flux is distributed acrass a gaussian around the + * center of an input pixel * * p: structure containing options, input, and output */ static int -do_kernel_gaussian(struct driz_param_t* p) { +do_kernel_gaussian(struct driz_param_t *p) { struct scanner s; integer_t bv, i, j, ii, jj, nxi, nxa, nyi, nya, nhit; integer_t osize[2]; float vc, d, dow; double gaussian_efac, gaussian_es; - double pfo, ac, scale2, xxi, xxa, yyi, yya, w, ddx, ddy, r2, dover; + double pfo, ac, scale2, xxi, xxa, yyi, yya, w, ddx, ddy, r2, dover; const double nsig = 2.5; int xmin, xmax, ymin, ymax, n; @@ -356,7 +367,7 @@ do_kernel_gaussian(struct driz_param_t* p) { scale2 = p->scale * p->scale; bv = compute_bit_value(p->uuid); - gaussian_efac = (2.3548*2.3548) * scale2 * ac / 2.0; + gaussian_efac = (2.3548 * 2.3548) * scale2 * ac / 2.0; gaussian_es = gaussian_efac / M_PI; if (init_image_scanner(p, &s, &ymin, &ymax)) return 1; @@ -376,8 +387,8 @@ do_kernel_gaussian(struct driz_param_t* p) { p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); break; } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) + // pixel centered on y is outside of scanner's limits or image [0, + // height - 1] OR: limits (x1, x2) are equal (line width is 0) p->nmiss += (p->xmax - p->xmin); ++p->nskip; continue; @@ -400,9 +411,9 @@ do_kernel_gaussian(struct driz_param_t* p) { yya = oy + pfo; nxi = MAX(fortran_round(xxi), 0); - nxa = MIN(fortran_round(xxa), osize[0]-1); + nxa = MIN(fortran_round(xxa), osize[0] - 1); nyi = MAX(fortran_round(yyi), 0); - nya = MIN(fortran_round(yya), osize[1]-1); + nya = MIN(fortran_round(yya), osize[1] - 1); nhit = 0; @@ -410,7 +421,8 @@ do_kernel_gaussian(struct driz_param_t* p) { d = get_pixel(p->data, i, j) * scale2; /* Scale the weighting mask by the scale factor and inversely by - the Jacobian to ensure conservation of weight in the output */ + the Jacobian to ensure conservation of weight in the output + */ if (p->weights) { w = get_pixel(p->weights, i, j) * p->weight_scale; } else { @@ -423,7 +435,7 @@ do_kernel_gaussian(struct driz_param_t* p) { for (ii = nxi; ii <= nxa; ++ii) { ddx = ox - (double)ii; /* Radial distance */ - r2 = ddx*ddx + ddy*ddy; + r2 = ddx * ddx + ddy * ddy; /* Weight is a scaled Gaussian function of radial distance */ @@ -435,8 +447,8 @@ do_kernel_gaussian(struct driz_param_t* p) { vc = get_pixel(p->output_counts, ii, jj); dow = (float)dover * w; - /* If we are create or modifying the context image, we do so - here. */ + /* If we are create or modifying the context image, we + do so here. */ if (p->output_context && dow > 0.0) { set_bit(p->output_context, ii, jj, bv); } @@ -449,21 +461,22 @@ do_kernel_gaussian(struct driz_param_t* p) { } /* Count cases where the pixel is off the output image */ - if (nhit == 0) ++ p->nmiss; + if (nhit == 0) ++p->nmiss; } } return 0; } -/** -------------------------------------------------------------------------------------------------- - * This kernel assumes flux of input pixel is distributed according to lanczos function +/** --------------------------------------------------------------------------- + * This kernel assumes flux of input pixel is distributed according to lanczos + * function * * p: structure containing options, input, and output */ static int -do_kernel_lanczos(struct driz_param_t* p) { +do_kernel_lanczos(struct driz_param_t *p) { struct scanner s; integer_t bv, i, j, ii, jj, nxi, nxa, nyi, nya, nhit, ix, iy; integer_t osize[2]; @@ -511,8 +524,8 @@ do_kernel_lanczos(struct driz_param_t* p) { p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); break; } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) + // pixel centered on y is outside of scanner's limits or image [0, + // height - 1] OR: limits (x1, x2) are equal (line width is 0) p->nmiss += (p->xmax - p->xmin); ++p->nskip; continue; @@ -532,9 +545,9 @@ do_kernel_lanczos(struct driz_param_t* p) { yya = yy - dy + pfo; nxi = MAX(fortran_round(xxi), 0); - nxa = MIN(fortran_round(xxa), osize[0]-1); + nxa = MIN(fortran_round(xxa), osize[0] - 1); nyi = MAX(fortran_round(yyi), 0); - nya = MIN(fortran_round(yya), osize[1]-1); + nya = MIN(fortran_round(yya), osize[1] - 1); nhit = 0; @@ -542,7 +555,8 @@ do_kernel_lanczos(struct driz_param_t* p) { d = get_pixel(p->data, i, j) * scale2; /* Scale the weighting mask by the scale factor and inversely by - the Jacobian to ensure conservation of weight in the output */ + the Jacobian to ensure conservation of weight in the output + */ if (p->weights) { w = get_pixel(p->weights, i, j) * p->weight_scale; } else { @@ -553,10 +567,15 @@ do_kernel_lanczos(struct driz_param_t* p) { for (jj = nyi; jj <= nya; ++jj) { for (ii = nxi; ii <= nxa; ++ii) { /* X and Y offsets */ - ix = fortran_round(fabs(xx - (double)ii) * lanczos.sdp) + 1; - iy = fortran_round(fabs(yy - (double)jj) * lanczos.sdp) + 1; - - /* Weight is product of Lanczos function values in X and Y */ + ix = + fortran_round(fabs(xx - (double)ii) * lanczos.sdp) + + 1; + iy = + fortran_round(fabs(yy - (double)jj) * lanczos.sdp) + + 1; + + /* Weight is product of Lanczos function values in X and + * Y */ dover = lanczos.lut[ix] * lanczos.lut[iy]; /* Count the hits */ @@ -565,8 +584,8 @@ do_kernel_lanczos(struct driz_param_t* p) { vc = get_pixel(p->output_counts, ii, jj); dow = (float)(dover * w); - /* If we are create or modifying the context image, we do so - here. */ + /* If we are create or modifying the context image, we + do so here. */ if (p->output_context && dow > 0.0) { set_bit(p->output_context, ii, jj, bv); } @@ -579,7 +598,7 @@ do_kernel_lanczos(struct driz_param_t* p) { } /* Count cases where the pixel is off the output image */ - if (nhit == 0) ++ p->nmiss; + if (nhit == 0) ++p->nmiss; } } @@ -589,15 +608,16 @@ do_kernel_lanczos(struct driz_param_t* p) { return 0; } -/** -------------------------------------------------------------------------------------------------- - * This kernel assumes the input flux is evenly distributed over a rectangle whose sides are - * aligned with the ouput pixel. Called turbo because it is fast, but approximate. +/** --------------------------------------------------------------------------- + * This kernel assumes the input flux is evenly distributed over a rectangle + * whose sides are aligned with the ouput pixel. Called turbo because it is + * fast, but approximate. * * p: structure containing options, input, and output */ static int -do_kernel_turbo(struct driz_param_t* p) { +do_kernel_turbo(struct driz_param_t *p) { struct scanner s; integer_t bv, i, j, ii, jj, nxi, nxa, nyi, nya, nhit, iis, iie, jjs, jje; integer_t osize[2]; @@ -630,8 +650,8 @@ do_kernel_turbo(struct driz_param_t* p) { p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); break; } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) + // pixel centered on y is outside of scanner's limits or image [0, + // height - 1] OR: limits (x1, x2) are equal (line width is 0) p->nmiss += (p->xmax - p->xmin); ++p->nskip; continue; @@ -657,10 +677,12 @@ do_kernel_turbo(struct driz_param_t* p) { nxa = fortran_round(xxa); nyi = fortran_round(yyi); nya = fortran_round(yya); - iis = MAX(nxi, 0); /* Needed to be set to 0 to avoid edge effects */ - iie = MIN(nxa, osize[0]-1); - jjs = MAX(nyi, 0); /* Needed to be set to 0 to avoid edge effects */ - jje = MIN(nya, osize[1]-1); + iis = MAX(nxi, + 0); /* Needed to be set to 0 to avoid edge effects */ + iie = MIN(nxa, osize[0] - 1); + jjs = MAX(nyi, + 0); /* Needed to be set to 0 to avoid edge effects */ + jje = MIN(nya, osize[1] - 1); nhit = 0; @@ -668,7 +690,8 @@ do_kernel_turbo(struct driz_param_t* p) { d = get_pixel(p->data, i, j) * (float)scale2; /* Scale the weighting mask by the scale factor and inversely by - the Jacobian to ensure conservation of weight in the output. */ + the Jacobian to ensure conservation of weight in the output. + */ if (p->weights) { w = get_pixel(p->weights, i, j) * p->weight_scale; } else { @@ -707,7 +730,7 @@ do_kernel_turbo(struct driz_param_t* p) { } /* Count cases where the pixel is off the output image */ - if (nhit == 0) ++ p->nmiss; + if (nhit == 0) ++p->nmiss; } } @@ -715,16 +738,17 @@ do_kernel_turbo(struct driz_param_t* p) { return 0; } -/** -------------------------------------------------------------------------------------------------- - * This module does the actual mapping of input flux to output images. It works by calculating the - * positions of the four corners of a quadrilateral on the output grid corresponding to the corners - * of the input pixel and then working out exactly how much of each pixel in the output is covered, or not. +/** --------------------------------------------------------------------------- + * This module does the actual mapping of input flux to output images. It works + * by calculating the positions of the four corners of a quadrilateral on the + * output grid corresponding to the corners of the input pixel and then working + * out exactly how much of each pixel in the output is covered, or not. * * p: structure containing options, input, and output */ int -do_kernel_square(struct driz_param_t* p) { +do_kernel_square(struct driz_param_t *p) { integer_t bv, i, j, ii, jj, min_ii, max_ii, min_jj, max_jj, nhit; integer_t osize[2]; float scale2, vc, d, dow; @@ -758,8 +782,8 @@ do_kernel_square(struct driz_param_t* p) { p->nmiss += (ymax + 1 - j) * (p->xmax - p->xmin); break; } else if (n == 2 || n == 3) { - // pixel centered on y is outside of scanner's limits or image [0, height - 1] - // OR: limits (x1, x2) are equal (line width is 0) + // pixel centered on y is outside of scanner's limits or image [0, + // height - 1] OR: limits (x1, x2) are equal (line width is 0) p->nmiss += (p->xmax - p->xmin); ++p->nskip; continue; @@ -770,17 +794,18 @@ do_kernel_square(struct driz_param_t* p) { /* Set the input corner positions */ - yin[1] = yin[0] = (double) j + dh; - yin[3] = yin[2] = (double) j - dh; + yin[1] = yin[0] = (double)j + dh; + yin[3] = yin[2] = (double)j - dh; for (i = xmin; i <= xmax; ++i) { nhit = 0; - xin[3] = xin[0] = (double) i - dh; - xin[2] = xin[1] = (double) i + dh; + xin[3] = xin[0] = (double)i - dh; + xin[2] = xin[1] = (double)i + dh; for (ii = 0; ii < 4; ++ii) { - if (interpolate_point(p, xin[ii], yin[ii], xout + ii, yout + ii)) { + if (interpolate_point(p, xin[ii], yin[ii], xout + ii, + yout + ii)) { goto _miss; } } @@ -795,8 +820,12 @@ do_kernel_square(struct driz_param_t* p) { if (jaco < 0.0) { jaco *= -1.0; /* Swap */ - tem = xout[1]; xout[1] = xout[3]; xout[3] = tem; - tem = yout[1]; yout[1] = yout[3]; yout[3] = tem; + tem = xout[1]; + xout[1] = xout[3]; + xout[3] = tem; + tem = yout[1]; + yout[1] = yout[3]; + yout[3] = tem; } /* Allow for stretching because of scale change */ @@ -812,9 +841,9 @@ do_kernel_square(struct driz_param_t* p) { /* Loop over output pixels which could be affected */ min_jj = MAX(fortran_round(min_doubles(yout, 4)), 0); - max_jj = MIN(fortran_round(max_doubles(yout, 4)), osize[1]-1); + max_jj = MIN(fortran_round(max_doubles(yout, 4)), osize[1] - 1); min_ii = MAX(fortran_round(min_doubles(xout, 4)), 0); - max_ii = MIN(fortran_round(max_doubles(xout, 4)), osize[0]-1); + max_ii = MIN(fortran_round(max_doubles(xout, 4)), osize[0] - 1); for (jj = min_jj; jj <= max_jj; ++jj) { for (ii = min_ii; ii <= max_ii; ++ii) { @@ -844,10 +873,10 @@ do_kernel_square(struct driz_param_t* p) { } } - /* Count cases where the pixel is off the output image */ - _miss: + /* Count cases where the pixel is off the output image */ + _miss: if (nhit == 0) { - ++ p->nmiss; + ++p->nmiss; } } } @@ -856,30 +885,25 @@ do_kernel_square(struct driz_param_t* p) { return 0; } -/** -------------------------------------------------------------------------------------------------- - * The user selects a kernel to use for drizzling from a function in the following tables - * The kernels differ in how the flux inside a single pixel is allocated: evenly spread - * across the pixel, concentrated at the central point, or by some other function. +/** --------------------------------------------------------------------------- + * The user selects a kernel to use for drizzling from a function in the + * following tables The kernels differ in how the flux inside a single pixel is + * allocated: evenly spread across the pixel, concentrated at the central point, + * or by some other function. */ -static kernel_handler_t -kernel_handler_map[] = { - do_kernel_square, - do_kernel_gaussian, - do_kernel_point, - do_kernel_turbo, - do_kernel_lanczos, - do_kernel_lanczos -}; - -/** -------------------------------------------------------------------------------------------------- +static kernel_handler_t kernel_handler_map[] = { + do_kernel_square, do_kernel_gaussian, do_kernel_point, + do_kernel_turbo, do_kernel_lanczos, do_kernel_lanczos}; + +/** --------------------------------------------------------------------------- * The executive function which calls the kernel which does the actual drizzling * * p: structure containing options, input, and output */ int -dobox(struct driz_param_t* p) { +dobox(struct driz_param_t *p) { kernel_handler_t kernel_handler = NULL; driz_log_message("starting dobox"); diff --git a/src/cdrizzlebox.h b/src/cdrizzlebox.h index 609e2799..374c3c78 100644 --- a/src/cdrizzlebox.h +++ b/src/cdrizzlebox.h @@ -19,18 +19,14 @@ In V1.6 this was simplified to use the DRIVAL routine and also to include some limited multi-kernel support. */ -integer_t -compute_bit_value(integer_t uuid); +integer_t compute_bit_value(integer_t uuid); -int -dobox(struct driz_param_t* p); +int dobox(struct driz_param_t *p); -double -compute_area(double is, double js, const double x[4], const double y[4]); +double compute_area(double is, double js, const double x[4], const double y[4]); -double -boxer(double is, double js, const double x[4], const double y[4]); +double boxer(double is, double js, const double x[4], const double y[4]); -typedef int (*kernel_handler_t)(struct driz_param_t*); +typedef int (*kernel_handler_t)(struct driz_param_t *); #endif /* CDRIZZLEBOX_H */ diff --git a/src/cdrizzlemap.h b/src/cdrizzlemap.h index 1024e30d..da9112e3 100644 --- a/src/cdrizzlemap.h +++ b/src/cdrizzlemap.h @@ -11,8 +11,8 @@ */ struct segment { - double point[2][2]; - int invalid; + double point[2][2]; + int invalid; }; // IMAGE_OUTLINE_NPTS - maximum number of vertices in the bounding polygon @@ -39,8 +39,8 @@ struct vertex { * */ struct polygon { - struct vertex v[2 * IMAGE_OUTLINE_NPTS]; /**< polygon vertices */ - int npv; /**< actual number of polygon vertices */ + struct vertex v[2 * IMAGE_OUTLINE_NPTS]; /**< polygon vertices */ + int npv; /**< actual number of polygon vertices */ }; /** edge structure. @@ -53,10 +53,11 @@ struct polygon { struct edge { struct vertex v1; /**< first vertex */ struct vertex v2; /**< second vertex */ - double m; /**< edge's slope */ - double b; /**< edge's interceipt */ - double c; /**< modified interceipt */ - int p; /**< edge's position: -1 for left-side edge and +1 for right-side edge */ + double m; /**< edge's slope */ + double b; /**< edge's interceipt */ + double c; /**< modified interceipt */ + int p; /**< edge's position: -1 for left-side edge and +1 for right-side + edge */ }; /** scanner structure. @@ -68,55 +69,55 @@ struct edge { * */ struct scanner { - struct edge left_edges[2 * IMAGE_OUTLINE_NPTS]; /**< left edges */ + struct edge left_edges[2 * IMAGE_OUTLINE_NPTS]; /**< left edges */ struct edge right_edges[2 * IMAGE_OUTLINE_NPTS]; /**< right edges */ - struct edge *left; /**< pointer to the current left edge; NULL when top polygon vertex was reached */ - struct edge *right; /**< pointer to the current right edge; NULL when top polygon vertex was reached */ - int nleft; /**< number of left edges */ - int nright; /**< number of right edges */ - double min_y; /**< minimum y-coordinate of all polygon vertices */ - double max_y; /**< maximum y-coordinate of all polygon vertices */ - int xmin; /**< min valid pixels' x-coord in pixmap (from bounding box carried over from driz_param_t) rounded to int */ - int xmax; /**< max valid pixels' x-coord in pixmap (from bounding box carried over from driz_param_t) rounded to int */ - int ymin; /**< min valid pixels' y-coord in pixmap (from bounding box carried over from driz_param_t) rounded to int */ - int ymax; /**< max valid pixels' y-coord in pixmap (from bounding box carried over from driz_param_t) rounded to int */ + struct edge *left; /**< pointer to the current left edge; NULL when top + polygon vertex was reached */ + struct edge *right; /**< pointer to the current right edge; NULL when top + polygon vertex was reached */ + int nleft; /**< number of left edges */ + int nright; /**< number of right edges */ + double min_y; /**< minimum y-coordinate of all polygon vertices */ + double max_y; /**< maximum y-coordinate of all polygon vertices */ + int xmin; /**< min valid pixels' x-coord in pixmap (from bounding box + carried over from driz_param_t) rounded to int */ + int xmax; /**< max valid pixels' x-coord in pixmap (from bounding box + carried over from driz_param_t) rounded to int */ + int ymin; /**< min valid pixels' y-coord in pixmap (from bounding box + carried over from driz_param_t) rounded to int */ + int ymax; /**< max valid pixels' y-coord in pixmap (from bounding box + carried over from driz_param_t) rounded to int */ // overlap_valid: 1 if polygon intersection and coord inversion worked; // 0 if computation of xmin, xmax, ymin, ymax has - // failed in which case they are carried over from driz_param_t. - int overlap_valid; /**< 1 if x/y min/max updated from polygon intersection and 0 if carried over from driz_param_t */ + // failed in which case they are carried over from + // driz_param_t. + int overlap_valid; /**< 1 if x/y min/max updated from polygon intersection + and 0 if carried over from driz_param_t */ }; -int -interpolate_point(struct driz_param_t *par, double xin, double yin, - double *xout, double *yout); +int interpolate_point(struct driz_param_t *par, double xin, double yin, + double *xout, double *yout); -int -map_point(struct driz_param_t *par, double xin, double yin, - double *xout, double *yout); +int map_point(struct driz_param_t *par, double xin, double yin, double *xout, + double *yout); -int -map_pixel(PyArrayObject *pixmap, int i, int j, double *x, double *y); +int map_pixel(PyArrayObject *pixmap, int i, int j, double *x, double *y); -int -shrink_image_section(PyArrayObject *pixmap, int *xmin, int *xmax, - int *ymin, int *ymax); +int shrink_image_section(PyArrayObject *pixmap, int *xmin, int *xmax, int *ymin, + int *ymax); -int -invert_pixmap(struct driz_param_t* par, double xout, double yout, - double *xin, double *yin); +int invert_pixmap(struct driz_param_t *par, double xout, double yout, + double *xin, double *yin); -int -clip_polygon_to_window(const struct polygon *p, const struct polygon *wnd, - struct polygon *cp); +int clip_polygon_to_window(const struct polygon *p, const struct polygon *wnd, + struct polygon *cp); -int -init_scanner(struct polygon *p, struct driz_param_t* par, struct scanner *s); +int init_scanner(struct polygon *p, struct driz_param_t *par, + struct scanner *s); -int -get_scanline_limits(struct scanner *s, int y, int *x1, int *x2); +int get_scanline_limits(struct scanner *s, int y, int *x1, int *x2); -int -init_image_scanner(struct driz_param_t* par, struct scanner *s, - int *ymin, int *ymax); +int init_image_scanner(struct driz_param_t *par, struct scanner *s, int *ymin, + int *ymax); #endif /* CDRIZZLEMAP_H */ diff --git a/src/cdrizzleutil.c b/src/cdrizzleutil.c index 4ea1f3cc..540ddb82 100644 --- a/src/cdrizzleutil.c +++ b/src/cdrizzleutil.c @@ -4,7 +4,7 @@ #include "cdrizzleutil.h" #include -#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ +#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ #include #include #include @@ -15,305 +15,282 @@ ERROR HANDLING */ void -driz_error_init(struct driz_error_t* error) { - assert(error); +driz_error_init(struct driz_error_t *error) { + assert(error); - error->last_message[0] = 0; + error->last_message[0] = 0; } int -driz_error_check(struct driz_error_t* error, const char* message, int test) { - if (! test) { - driz_error_set_message(error, message); - return 1; - } +driz_error_check(struct driz_error_t *error, const char *message, int test) { + if (!test) { + driz_error_set_message(error, message); + return 1; + } - return 0; + return 0; } void -driz_error_set_message(struct driz_error_t* error, const char* message) { - assert(error); - assert(message); +driz_error_set_message(struct driz_error_t *error, const char *message) { + assert(error); + assert(message); - strncpy(error->last_message, message, MAX_DRIZ_ERROR_LEN); + strncpy(error->last_message, message, MAX_DRIZ_ERROR_LEN); } void -driz_error_format_message(struct driz_error_t* error, const char* format, ...) { - /* See http://c-faq.com/varargs/vprintf.html - for an explanation of how all this variable length argument list stuff - works. */ - va_list argp; - - assert(error); - assert(format); - - va_start(argp, format); - (void)vsnprintf(error->last_message, MAX_DRIZ_ERROR_LEN, format, argp); - va_end(argp); +driz_error_format_message(struct driz_error_t *error, const char *format, ...) { + /* See http://c-faq.com/varargs/vprintf.html + for an explanation of how all this variable length argument list stuff + works. */ + va_list argp; + + assert(error); + assert(format); + + va_start(argp, format); + (void)vsnprintf(error->last_message, MAX_DRIZ_ERROR_LEN, format, argp); + va_end(argp); } -const char* -driz_error_get_message(struct driz_error_t* error) { - assert(error); +const char * +driz_error_get_message(struct driz_error_t *error) { + assert(error); - return error->last_message; + return error->last_message; } int -driz_error_is_set(struct driz_error_t* error) { - assert(error); +driz_error_is_set(struct driz_error_t *error) { + assert(error); - return error->last_message[0] != 0; + return error->last_message[0] != 0; } void -driz_error_unset(struct driz_error_t* error) { - assert(error); +driz_error_unset(struct driz_error_t *error) { + assert(error); - driz_error_init(error); + driz_error_init(error); } /***************************************************************** DATA TYPES */ void -driz_param_dump(struct driz_param_t* p) { - assert(p); - - printf("DRIZZLING PARAMETERS:\n" - "kernel: %s\n" - "pixel_fraction: %f\n" - "exposure_time: %f\n" - "weight_scale: %f\n" - "fill_value: %f\n" - "do_fill: %s\n" - "in_units: %s\n" - "out_units: %s\n" - "scale: %f\n", - kernel_enum2str(p->kernel), - p->pixel_fraction, - p->exposure_time, - p->weight_scale, - p->fill_value, - bool2str(p->do_fill), - unit_enum2str(p->in_units), - unit_enum2str(p->out_units), - p->scale - ); +driz_param_dump(struct driz_param_t *p) { + assert(p); + + printf( + "DRIZZLING PARAMETERS:\n" + "kernel: %s\n" + "pixel_fraction: %f\n" + "exposure_time: %f\n" + "weight_scale: %f\n" + "fill_value: %f\n" + "do_fill: %s\n" + "in_units: %s\n" + "out_units: %s\n" + "scale: %f\n", + kernel_enum2str(p->kernel), p->pixel_fraction, p->exposure_time, + p->weight_scale, p->fill_value, bool2str(p->do_fill), + unit_enum2str(p->in_units), unit_enum2str(p->out_units), p->scale); } void -driz_param_init(struct driz_param_t* p) { - assert(p); +driz_param_init(struct driz_param_t *p) { + assert(p); - /* Kernel shape and size */ - p->kernel = kernel_square; - p->pixel_fraction = 1.0; + /* Kernel shape and size */ + p->kernel = kernel_square; + p->pixel_fraction = 1.0; - /* Exposure time */ - p->exposure_time = 1.0; + /* Exposure time */ + p->exposure_time = 1.0; - /* Weight scale */ - p->weight_scale = 1.0; + /* Weight scale */ + p->weight_scale = 1.0; - /* Filling */ - p->fill_value = 0.0; - p->do_fill = 0; + /* Filling */ + p->fill_value = 0.0; + p->do_fill = 0; - /* CPS / Counts */ - p->in_units = unit_counts; - p->out_units = unit_counts; + /* CPS / Counts */ + p->in_units = unit_counts; + p->out_units = unit_counts; - p->scale = 1.0; + p->scale = 1.0; - /* Input data */ - p->data = NULL; - p->weights = NULL; - p->pixmap = NULL; + /* Input data */ + p->data = NULL; + p->weights = NULL; + p->pixmap = NULL; - /* Output data */ - p->output_data = NULL; - p->output_counts = NULL; - p->output_context = NULL; + /* Output data */ + p->output_data = NULL; + p->output_counts = NULL; + p->output_context = NULL; - p->nmiss = 0; - p->nskip = 0; - p->error = NULL; + p->nmiss = 0; + p->nskip = 0; + p->error = NULL; } /***************************************************************** STRING TO ENUMERATION CONVERSIONS */ -static const char* kernel_string_table[] = { - "square", - "gaussian", - "point", - "turbo", - "lanczos2", - "lanczos3", - NULL -}; - -static const char* unit_string_table[] = { - "counts", - "cps", - NULL -}; - -static const char* interp_string_table[] = { - "nearest", - "linear", - "poly3", - "poly5", - "spline3", - "sinc", - "lsinc", - "lan3", - "lan5", - NULL -}; - -static const char* bool_string_table[] = { - "FALSE", - "TRUE", - NULL -}; +static const char *kernel_string_table[] = { + "square", "gaussian", "point", "turbo", "lanczos2", "lanczos3", NULL}; + +static const char *unit_string_table[] = {"counts", "cps", NULL}; + +static const char *interp_string_table[] = { + "nearest", "linear", "poly3", "poly5", "spline3", + "sinc", "lsinc", "lan3", "lan5", NULL}; + +static const char *bool_string_table[] = {"FALSE", "TRUE", NULL}; static int -str2enum(const char* s, const char* table[], int* result, struct driz_error_t* error) { - const char** it = table; - - assert(s); - assert(table); - assert(result); - assert(error); - - while (*it != NULL) { - if (strncmp(s, *it, 32) == 0) { - *result = it - table; - return 0; +str2enum(const char *s, const char *table[], int *result, + struct driz_error_t *error) { + const char **it = table; + + assert(s); + assert(table); + assert(result); + assert(error); + + while (*it != NULL) { + if (strncmp(s, *it, 32) == 0) { + *result = it - table; + return 0; + } + ++it; } - ++it; - } - return 1; + return 1; } int -kernel_str2enum(const char* s, enum e_kernel_t* result, struct driz_error_t* error) { - if (str2enum(s, kernel_string_table, (int *)result, error)) { - driz_error_format_message(error, "Unknown kernel type '%s'", s); - return 1; - } +kernel_str2enum(const char *s, enum e_kernel_t *result, + struct driz_error_t *error) { + if (str2enum(s, kernel_string_table, (int *)result, error)) { + driz_error_format_message(error, "Unknown kernel type '%s'", s); + return 1; + } - return 0; + return 0; } int -unit_str2enum(const char* s, enum e_unit_t* result, struct driz_error_t* error) { - if (str2enum(s, unit_string_table, (int *)result, error)) { - driz_error_format_message(error, "Unknown unit type '%s'", s); - return 1; - } +unit_str2enum(const char *s, enum e_unit_t *result, + struct driz_error_t *error) { + if (str2enum(s, unit_string_table, (int *)result, error)) { + driz_error_format_message(error, "Unknown unit type '%s'", s); + return 1; + } - return 0; + return 0; } int -interp_str2enum(const char* s, enum e_interp_t* result, struct driz_error_t* error) { - if (str2enum(s, interp_string_table, (int *)result, error)) { - driz_error_format_message(error, "Unknown interp type '%s'", s); - return 1; - } +interp_str2enum(const char *s, enum e_interp_t *result, + struct driz_error_t *error) { + if (str2enum(s, interp_string_table, (int *)result, error)) { + driz_error_format_message(error, "Unknown interp type '%s'", s); + return 1; + } - return 0; + return 0; } -const char* +const char * kernel_enum2str(enum e_kernel_t value) { - assert(value >= 0 && value < kernel_LAST); + assert(value >= 0 && value < kernel_LAST); - return kernel_string_table[value]; + return kernel_string_table[value]; } -const char* +const char * unit_enum2str(enum e_unit_t value) { - assert(value >= 0 && value < 2); + assert(value >= 0 && value < 2); - return unit_string_table[value]; + return unit_string_table[value]; } -const char* +const char * interp_enum2str(enum e_interp_t value) { - assert(value >= 0 && value < interp_LAST); + assert(value >= 0 && value < interp_LAST); - return interp_string_table[value]; + return interp_string_table[value]; } -const char* +const char * bool2str(bool_t value) { - return bool_string_table[value ? 1 : 0]; + return bool_string_table[value ? 1 : 0]; } /***************************************************************** NUMERICAL UTILITIES */ void -create_lanczos_lut(const int kernel_order, const size_t npix, - const float del, float* lanczos_lut) { - integer_t i; - const float forder = (float)kernel_order; - float poff; - - assert(lanczos_lut); - assert(kernel_order < 6); - - /* Set the first value to avoid arithmetic problems */ - lanczos_lut[0] = 1.0; - - for (i = 1; i < npix; ++i) { - poff = M_PI * (float)i * del; - if (poff < M_PI * forder) { - lanczos_lut[i] = sin(poff) / poff * sin(poff / forder) / (poff / forder); - } else { - lanczos_lut[i] = 0.0; +create_lanczos_lut(const int kernel_order, const size_t npix, const float del, + float *lanczos_lut) { + integer_t i; + const float forder = (float)kernel_order; + float poff; + + assert(lanczos_lut); + assert(kernel_order < 6); + + /* Set the first value to avoid arithmetic problems */ + lanczos_lut[0] = 1.0; + + for (i = 1; i < npix; ++i) { + poff = M_PI * (float)i * del; + if (poff < M_PI * forder) { + lanczos_lut[i] = + sin(poff) / poff * sin(poff / forder) / (poff / forder); + } else { + lanczos_lut[i] = 0.0; + } } - } } void -put_fill(struct driz_param_t* p, const float fill_value) { - integer_t i, j, osize[2]; - - assert(p); - get_dimensions(p->output_data, osize); - for (j = 0; j < osize[1]; ++j) { - for (i = 0; i < osize[0]; ++i) { - if (oob_pixel(p->output_counts, i, j)) { - driz_error_format_message(p->error, "OOB in output_counts[%d,%d]", i, j); - return; - - } else if (oob_pixel(p->output_data, i, j)) { - driz_error_format_message(p->error, "OOB in output_data[%d,%d]", i, j); - return; - - } else if (get_pixel(p->output_counts, i, j) == 0.0) { - set_pixel(p->output_data, i, j, fill_value); - } +put_fill(struct driz_param_t *p, const float fill_value) { + integer_t i, j, osize[2]; + + assert(p); + get_dimensions(p->output_data, osize); + for (j = 0; j < osize[1]; ++j) { + for (i = 0; i < osize[0]; ++i) { + if (oob_pixel(p->output_counts, i, j)) { + driz_error_format_message(p->error, + "OOB in output_counts[%d,%d]", i, j); + return; + + } else if (oob_pixel(p->output_data, i, j)) { + driz_error_format_message(p->error, "OOB in output_data[%d,%d]", + i, j); + return; + + } else if (get_pixel(p->output_counts, i, j) == 0.0) { + set_pixel(p->output_data, i, j, fill_value); + } + } } - } } double mgf2(double lambda) { - double sig, sig2; + double sig, sig2; - sig = 1.0e7 / lambda; - sig2 = sig*sig; + sig = 1.0e7 / lambda; + sig2 = sig * sig; - return sqrt(1.0 + 2.590355e10/(5.312993e10-sig2) + - 4.4543708e9/(11.17083e9-sig2) + - 4.0838897e5/(1.766361e5-sig2)); + return sqrt(1.0 + 2.590355e10 / (5.312993e10 - sig2) + + 4.4543708e9 / (11.17083e9 - sig2) + + 4.0838897e5 / (1.766361e5 - sig2)); } diff --git a/src/cdrizzleutil.h b/src/cdrizzleutil.h index d230bc5f..04d48236 100644 --- a/src/cdrizzleutil.h +++ b/src/cdrizzleutil.h @@ -10,7 +10,7 @@ #include #include -#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ +#define _USE_MATH_DEFINES /* needed for MS Windows to define M_PI */ #include #if __STDC_VERSION__ >= 199901L #include @@ -23,31 +23,33 @@ #define MAX_DRIZ_ERROR_LEN 512 struct driz_error_t { - char last_message[MAX_DRIZ_ERROR_LEN]; + char last_message[MAX_DRIZ_ERROR_LEN]; }; -void driz_error_init(struct driz_error_t* error); -int driz_error_check(struct driz_error_t* error, const char* message, int test); -void driz_error_set_message(struct driz_error_t* error, const char* message); -void driz_error_format_message(struct driz_error_t* error, const char* format, ...); -const char* driz_error_get_message(struct driz_error_t* error); -int driz_error_is_set(struct driz_error_t* error); -void driz_error_unset(struct driz_error_t* error); +void driz_error_init(struct driz_error_t *error); +int driz_error_check(struct driz_error_t *error, const char *message, int test); +void driz_error_set_message(struct driz_error_t *error, const char *message); +void driz_error_format_message(struct driz_error_t *error, const char *format, + ...); +const char *driz_error_get_message(struct driz_error_t *error); +int driz_error_is_set(struct driz_error_t *error); +void driz_error_unset(struct driz_error_t *error); /***************************************************************** CONVENIENCE MACROS */ #if !defined(MIN) - #define MIN(a, b) (((a) < (b)) ? (a) : (b)) +#define MIN(a, b) (((a) < (b)) ? (a) : (b)) #endif #if !defined(MAX) - #define MAX(a, b) (((a) > (b)) ? (a) : (b)) +#define MAX(a, b) (((a) > (b)) ? (a) : (b)) #endif -#define CLAMP(x, low, high) (((x) > (high)) ? (high) : (((x) < (low)) ? (low) : (x))) +#define CLAMP(x, low, high) \ + (((x) > (high)) ? (high) : (((x) < (low)) ? (low) : (x))) #define CLAMP_ABOVE(x, low) (((x) < low) ? (low) : (x)) -#define CLAMP_BELOW(x, high) (((x) > high) ? (high) : (x)) +#define CLAMP_BELOW(x, high) (((x) > high) ? (high) : (x)) #ifdef __GNUC__ #define UNUSED_PARAM __attribute__((unused)) @@ -58,7 +60,7 @@ void driz_error_unset(struct driz_error_t* error); #define MAX_DOUBLE 1.7976931348623158e+308 #define MIN_DOUBLE 2.2250738585072014e-308 -#define MAX_COEFFS 128 +#define MAX_COEFFS 128 #define COEFF_OFFSET 100 #undef TRUE @@ -78,86 +80,82 @@ typedef unsigned char bool_t; #endif enum e_kernel_t { - kernel_square, - kernel_gaussian, - kernel_point, - kernel_turbo, - kernel_lanczos2, - kernel_lanczos3, - kernel_LAST + kernel_square, + kernel_gaussian, + kernel_point, + kernel_turbo, + kernel_lanczos2, + kernel_lanczos3, + kernel_LAST }; -enum e_unit_t { - unit_counts, - unit_cps -}; +enum e_unit_t { unit_counts, unit_cps }; enum e_interp_t { - interp_nearest, - interp_bilinear, - interp_poly3, - interp_poly5, - interp_spline3, - interp_sinc, - interp_lsinc, - interp_lanczos3, - interp_lanczos5, - interp_LAST + interp_nearest, + interp_bilinear, + interp_poly3, + interp_poly5, + interp_spline3, + interp_sinc, + interp_lsinc, + interp_lanczos3, + interp_lanczos5, + interp_LAST }; /* Lanczos values */ struct lanczos_param_t { - size_t nlut; - float* lut; - double sdp; - integer_t nbox; - float space; - float misval; + size_t nlut; + float *lut; + double sdp; + integer_t nbox; + float space; + float misval; }; struct driz_param_t { - /* Options */ - enum e_kernel_t kernel; /* Kernel shape and size */ - double pixel_fraction; /* was: PIXFRAC */ - float exposure_time; /* Exposure time was: EXPIN */ - float weight_scale; /* Weight scale was: WTSCL */ - float fill_value; /* Filling was: FILVAL */ - bool_t do_fill; /* was: FILL */ - enum e_unit_t in_units; /* CPS / counts was: INCPS, either counts or CPS */ - enum e_unit_t out_units; /* CPS / counts was: INCPS, either counts or CPS */ - integer_t uuid; /* was: UNIQID */ - - /* Scaling */ - double scale; - - /* Image subset */ - integer_t xmin; - integer_t xmax; - integer_t ymin; - integer_t ymax; - - /* Blotting-specific parameters */ - enum e_interp_t interpolation; /* was INTERP */ - float ef; - float misval; - float sinscl; - float kscale; - - /* Input images */ - PyArrayObject *data; - PyArrayObject *weights; - PyArrayObject *pixmap; - - /* Output images */ - PyArrayObject *output_data; - PyArrayObject *output_counts; /* was: COU */ - PyArrayObject *output_context; /* was: CONTIM */ - - /* Other output */ - integer_t nmiss; - integer_t nskip; - struct driz_error_t* error; - + /* Options */ + enum e_kernel_t kernel; /* Kernel shape and size */ + double pixel_fraction; /* was: PIXFRAC */ + float exposure_time; /* Exposure time was: EXPIN */ + float weight_scale; /* Weight scale was: WTSCL */ + float fill_value; /* Filling was: FILVAL */ + bool_t do_fill; /* was: FILL */ + enum e_unit_t in_units; /* CPS / counts was: INCPS, either counts or CPS */ + enum e_unit_t out_units; /* CPS / counts was: INCPS, either counts or CPS */ + integer_t uuid; /* was: UNIQID */ + + /* Scaling */ + double scale; + + /* Image subset */ + integer_t xmin; + integer_t xmax; + integer_t ymin; + integer_t ymax; + + /* Blotting-specific parameters */ + enum e_interp_t interpolation; /* was INTERP */ + float ef; + float misval; + float sinscl; + float kscale; + + /* Input images */ + PyArrayObject *data; + PyArrayObject *weights; + PyArrayObject *pixmap; + + /* Output images */ + PyArrayObject *output_data; + PyArrayObject *output_counts; /* was: COU */ + PyArrayObject *output_context; /* was: CONTIM */ + + /* Other output */ + integer_t nmiss; + integer_t nskip; + struct driz_error_t *error; }; /** @@ -166,12 +164,9 @@ values, mostly zeroes. Note, these are not *meaningful* values, just ones that help with memory management etc. It is up to users of the struct, e.g. cdrizzle_, to fill the struct with valid parameters. */ -void -driz_param_init(struct driz_param_t* p); - -void -driz_param_dump(struct driz_param_t* p); +void driz_param_init(struct driz_param_t *p); +void driz_param_dump(struct driz_param_t *p); /****************************************************************************/ /* LOGGING */ @@ -179,10 +174,9 @@ driz_param_dump(struct driz_param_t* p); #ifdef LOGGING extern FILE *driz_log_handle; - -static inline_macro FILE* +static inline_macro FILE * driz_log_init(FILE *handle) { - const char* dirs[] = {"TMPDIR", "TMP", "TEMP", "TEMPDIR"}; + const char *dirs[] = {"TMPDIR", "TMP", "TEMP", "TEMPDIR"}; int i; char *p; char buf[1024]; @@ -210,7 +204,6 @@ driz_log_init(FILE *handle) { return handle; } - static inline_macro int driz_log_close(FILE *handle) { if (handle) { @@ -219,7 +212,7 @@ driz_log_close(FILE *handle) { } static inline_macro int -driz_log_message(const char* message) { +driz_log_message(const char *message) { if (!driz_log_handle) { driz_log_handle = driz_log_init(driz_log_handle); if (!driz_log_handle) { @@ -238,8 +231,8 @@ driz_log_idem(void *ptr) { return ptr; } -#define driz_log_init(handle) driz_log_idem(handle) -#define driz_log_close(handle) driz_log_idem(handle) +#define driz_log_init(handle) driz_log_idem(handle) +#define driz_log_close(handle) driz_log_idem(handle) #define driz_log_message(message) driz_log_idem(message) #endif @@ -251,19 +244,18 @@ driz_log_idem(void *ptr) { static inline_macro void get_dimensions(PyArrayObject *image, integer_t size[2]) { + npy_intp *ndim = PyArray_DIMS(image); - npy_intp *ndim = PyArray_DIMS(image); - - /* Put dimensions in xy order */ - size[0] = ndim[1]; - size[1] = ndim[0]; + /* Put dimensions in xy order */ + size[0] = ndim[1]; + size[1] = ndim[0]; - return; + return; } -static inline_macro double* +static inline_macro double * get_pixmap(PyArrayObject *pixmap, integer_t xpix, integer_t ypix) { - return (double*) PyArray_GETPTR3(pixmap, ypix, xpix, 0); + return (double *)PyArray_GETPTR3(pixmap, ypix, xpix, 0); } #if defined(LOGGING) && defined(CHECK_OOB) @@ -273,8 +265,8 @@ oob_pixel(PyArrayObject *image, integer_t xpix, integer_t ypix) { char buffer[64]; npy_intp *ndim = PyArray_DIMS(image); if ((xpix < 0 || xpix >= ndim[1]) || (ypix < 0 || ypix >= ndim[0])) { - sprintf(buffer, "Point [%d,%d] is outside of [%d, %d]", - xpix, ypix, (int) ndim[1], (int) ndim[0]); + sprintf(buffer, "Point [%d,%d] is outside of [%d, %d]", xpix, ypix, + (int)ndim[1], (int)ndim[0]); driz_log_message(buffer); return 1; } @@ -290,64 +282,62 @@ oob_pixel(PyArrayObject *image, integer_t xpix, integer_t ypix) { static inline_macro float get_pixel(PyArrayObject *image, integer_t xpix, integer_t ypix) { - return *(float*) PyArray_GETPTR2(image, ypix, xpix); + return *(float *)PyArray_GETPTR2(image, ypix, xpix); } static inline_macro float get_pixel_at_pos(PyArrayObject *image, integer_t pos) { - float *imptr; - imptr = (float *) PyArray_DATA(image); - return imptr[pos]; + float *imptr; + imptr = (float *)PyArray_DATA(image); + return imptr[pos]; } static inline_macro void set_pixel(PyArrayObject *image, integer_t xpix, integer_t ypix, double value) { - *(float*) PyArray_GETPTR2(image, ypix, xpix) = value; - return; + *(float *)PyArray_GETPTR2(image, ypix, xpix) = value; + return; } static inline_macro int -get_bit(PyArrayObject *image, integer_t xpix, integer_t ypix, integer_t bitval) { - integer_t value; - value = *(integer_t*) PyArray_GETPTR2(image, ypix, xpix) & bitval; - return value? 1 : 0; +get_bit(PyArrayObject *image, integer_t xpix, integer_t ypix, + integer_t bitval) { + integer_t value; + value = *(integer_t *)PyArray_GETPTR2(image, ypix, xpix) & bitval; + return value ? 1 : 0; } static inline_macro void -set_bit(PyArrayObject *image, integer_t xpix, integer_t ypix, integer_t bitval) { - *(integer_t*) PyArray_GETPTR2(image, ypix, xpix) |= bitval; - return; +set_bit(PyArrayObject *image, integer_t xpix, integer_t ypix, + integer_t bitval) { + *(integer_t *)PyArray_GETPTR2(image, ypix, xpix) |= bitval; + return; } static inline_macro void unset_bit(PyArrayObject *image, integer_t xpix, integer_t ypix) { - *(integer_t*) PyArray_GETPTR2(image, ypix, xpix) = 0; - return; + *(integer_t *)PyArray_GETPTR2(image, ypix, xpix) = 0; + return; } /***************************************************************** STRING TO ENUMERATION CONVERSIONS */ -int -kernel_str2enum(const char* s, enum e_kernel_t* result, struct driz_error_t* error); +int kernel_str2enum(const char *s, enum e_kernel_t *result, + struct driz_error_t *error); -int -unit_str2enum(const char* s, enum e_unit_t* result, struct driz_error_t* error); +int unit_str2enum(const char *s, enum e_unit_t *result, + struct driz_error_t *error); -int -interp_str2enum(const char* s, enum e_interp_t* result, struct driz_error_t* error); +int interp_str2enum(const char *s, enum e_interp_t *result, + struct driz_error_t *error); -const char* -kernel_enum2str(enum e_kernel_t value); +const char *kernel_enum2str(enum e_kernel_t value); -const char* -unit_enum2str(enum e_unit_t value); +const char *unit_enum2str(enum e_unit_t value); -const char* -interp_enum2str(enum e_interp_t value); +const char *interp_enum2str(enum e_interp_t value); -const char* -bool2str(bool_t value); +const char *bool2str(bool_t value); /***************************************************************** NUMERICAL UTILITIES @@ -366,19 +356,16 @@ Note that no checking is done to see whether the values are sensible. was: FILALU */ -void -create_lanczos_lut(const int kernel_order, const size_t npix, - const float del, float* lanczos_lut); +void create_lanczos_lut(const int kernel_order, const size_t npix, + const float del, float *lanczos_lut); -void -put_fill(struct driz_param_t* p, const float fill_value); +void put_fill(struct driz_param_t *p, const float fill_value); /** Calculate the refractive index of MgF2 for a given C wavelength (in nm) using the formula given by Trauger (1995) */ -double -mgf2(double lambda); +double mgf2(double lambda); /** Weighted sum of 2 real vectors. @@ -386,19 +373,17 @@ Weighted sum of 2 real vectors. was: WSUMR */ static inline_macro void -weighted_sum_vectors(const integer_t npix, - const float* a /*[npix]*/, const float w1, - const float* b /*[npix]*/, const float w2, +weighted_sum_vectors(const integer_t npix, const float *a /*[npix]*/, + const float w1, const float *b /*[npix]*/, const float w2, /* Output arguments */ - float* c /*[npix]*/) { - float* c_end = c + npix; + float *c /*[npix]*/) { + float *c_end = c + npix; - assert(a); - assert(b); - assert(c); + assert(a); + assert(b); + assert(c); - while(c != c_end) - *(c++) = *(a++) * w1 + *(b++) * w2; + while (c != c_end) *(c++) = *(a++) * w1 + *(b++) * w2; } /** @@ -406,27 +391,25 @@ weighted_sum_vectors(const integer_t npix, */ static inline_macro integer_t fortran_round(const double x) { - return (x >= 0) ? (integer_t)floor(x + .5) : (integer_t)-floor(.5 - x); + return (x >= 0) ? (integer_t)floor(x + .5) : (integer_t)-floor(.5 - x); } static inline_macro double -min_doubles(const double* a, const integer_t size) { - const double* end = a + size; - double value = MAX_DOUBLE; - for ( ; a != end; ++a) - if (*a < value) - value = *a; - return value; +min_doubles(const double *a, const integer_t size) { + const double *end = a + size; + double value = MAX_DOUBLE; + for (; a != end; ++a) + if (*a < value) value = *a; + return value; } static inline_macro double -max_doubles(const double* a, const integer_t size) { - const double* end = a + size; - double value = MIN_DOUBLE; - for ( ; a != end; ++a) - if (*a > value) - value = *a; - return value; +max_doubles(const double *a, const integer_t size) { + const double *end = a + size; + double value = MIN_DOUBLE; + for (; a != end; ++a) + if (*a > value) value = *a; + return value; } /** @@ -445,20 +428,20 @@ as this is physically meaningless. @param[out] yo The distorted y coordinate */ static inline_macro void -rad3(const double x, const double y, const double* co, +rad3(const double x, const double y, const double *co, /* Output parameters */ - double* xo, double* yo) { - double r, f; + double *xo, double *yo) { + double r, f; - assert(co); - assert(xo); - assert(yo); + assert(co); + assert(xo); + assert(yo); - r = sqrt(x*x + y*y); + r = sqrt(x * x + y * y); - f = 1.0 + co[0] + co[1]*r + co[2]*r*r; - *xo = f*x; - *yo = f*y; + f = 1.0 + co[0] + co[1] * r + co[2] * r * r; + *xo = f * x; + *yo = f * y; } #endif /* CDRIZZLEUTIL_H */ diff --git a/src/driz_portability.h b/src/driz_portability.h index 07292766..8faa8290 100644 --- a/src/driz_portability.h +++ b/src/driz_portability.h @@ -2,8 +2,8 @@ #define inline_macro __inline #else /* -* assume gcc for now -*/ + * assume gcc for now + */ #define inline_macro inline #endif diff --git a/src/tests/.clang-format b/src/tests/.clang-format new file mode 100644 index 00000000..2225c245 --- /dev/null +++ b/src/tests/.clang-format @@ -0,0 +1,5 @@ +--- +DisableFormat: true +SortIncludes: Never + +... diff --git a/src/tests/drizzletest.h b/src/tests/drizzletest.h index 6cccede6..ec07db72 100644 --- a/src/tests/drizzletest.h +++ b/src/tests/drizzletest.h @@ -6,9 +6,10 @@ #endif #include -int do_kernel_square(struct driz_param_t* p); +int do_kernel_square(struct driz_param_t *p); void set_test_arrays(PyArrayObject *dat, PyArrayObject *wei, PyArrayObject *map, - PyArrayObject *odat, PyArrayObject *ocnt, PyArrayObject *ocon); + PyArrayObject *odat, PyArrayObject *ocnt, + PyArrayObject *ocon); -int utest_cdrizzle(int argc, char* argv[]); +int utest_cdrizzle(int argc, char *argv[]); diff --git a/src/tests/utest_cdrizzle.c b/src/tests/utest_cdrizzle.c index b32e2fd9..c7971c69 100644 --- a/src/tests/utest_cdrizzle.c +++ b/src/tests/utest_cdrizzle.c @@ -21,7 +21,6 @@ #include "cdrizzleutil.h" #include "drizzletest.h" - FILE *logptr = NULL; static integer_t image_size[2]; static PyArrayObject *test_data; @@ -34,13 +33,8 @@ static PyArrayObject *test_context; static char log_file[] = ""; void -set_test_arrays(PyArrayObject *dat, - PyArrayObject *wei, - PyArrayObject *map, - PyArrayObject *odat, - PyArrayObject *ocnt, - PyArrayObject *ocon) { - +set_test_arrays(PyArrayObject *dat, PyArrayObject *wei, PyArrayObject *map, + PyArrayObject *odat, PyArrayObject *ocnt, PyArrayObject *ocon) { test_data = dat; test_weights = wei; test_pixmap = map; @@ -59,13 +53,13 @@ set_pixmap(struct driz_param_t *p, int xmin, int xmax, int ymin, int ymax) { ypix = ymin; for (j = ymin; j < ymax; j++) { - xpix = xmin; - for (i = xmin; i < xmax; i++) { + xpix = xmin; + for (i = xmin; i < xmax; i++) { get_pixmap(p->pixmap, i, j)[0] = xpix; get_pixmap(p->pixmap, i, j)[1] = ypix; xpix += 1.0; - } - ypix += 1.0; + } + ypix += 1.0; } return; @@ -79,19 +73,18 @@ init_pixmap(struct driz_param_t *p) { void stretch_pixmap(struct driz_param_t *p, double stretch) { - int i, j; double xpix, ypix; ypix = 0.0; for (j = 0; j < image_size[1]; j++) { - xpix = 0.0; - for (i= 0; i < image_size[0]; i++) { + xpix = 0.0; + for (i = 0; i < image_size[0]; i++) { get_pixmap(p->pixmap, i, j)[0] = xpix; get_pixmap(p->pixmap, i, j)[1] = stretch * ypix; xpix += 1.0; - } - ypix += 1.0; + } + ypix += 1.0; } return; @@ -102,10 +95,10 @@ nan_pixmap(struct driz_param_t *p) { int i, j; for (j = 0; j < image_size[1]; j++) { - for (i= 0; i < image_size[0]; i++) { + for (i = 0; i < image_size[0]; i++) { get_pixmap(p->pixmap, i, j)[0] = NPY_NAN; get_pixmap(p->pixmap, i, j)[1] = NPY_NAN; - } + } } return; @@ -115,25 +108,24 @@ void nan_pixel(struct driz_param_t *p, int xpix, int ypix) { int idim; for (idim = 0; idim < 2; ++idim) { - get_pixmap(p->pixmap, xpix, ypix)[idim] = NPY_NAN; + get_pixmap(p->pixmap, xpix, ypix)[idim] = NPY_NAN; } } void offset_pixmap(struct driz_param_t *p, double x_offset, double y_offset) { - int i, j; double xpix, ypix; ypix = y_offset; for (j = 0; j < image_size[1]; j++) { - xpix = x_offset; - for (i = 0; i < image_size[0]; i++) { + xpix = x_offset; + for (i = 0; i < image_size[0]; i++) { get_pixmap(p->pixmap, i, j)[0] = xpix; get_pixmap(p->pixmap, i, j)[1] = ypix; xpix += 1.0; - } - ypix += 1.0; + } + ypix += 1.0; } return; @@ -141,7 +133,7 @@ offset_pixmap(struct driz_param_t *p, double x_offset, double y_offset) { void fill_image(PyArrayObject *image, double value) { - npy_intp *ndim = PyArray_DIMS(image); + npy_intp *ndim = PyArray_DIMS(image); int ypix, xpix; for (ypix = 0; ypix < ndim[0]; ++ypix) { @@ -154,7 +146,7 @@ fill_image(PyArrayObject *image, double value) { } void -fill_image_block(PyArrayObject* image, double value, int lo, int hi) { +fill_image_block(PyArrayObject *image, double value, int lo, int hi) { int ypix, xpix; for (ypix = lo; ypix < hi; ++ypix) { @@ -167,7 +159,7 @@ fill_image_block(PyArrayObject* image, double value, int lo, int hi) { } void unset_context(PyArrayObject *context) { - npy_intp *ndim = PyArray_DIMS(context); + npy_intp *ndim = PyArray_DIMS(context); int ypix, xpix; for (ypix = 0; ypix < ndim[0]; ++ypix) { @@ -180,7 +172,7 @@ unset_context(PyArrayObject *context) { } void -print_image(char *title, PyArrayObject* image, int lo, int hi) { +print_image(char *title, PyArrayObject *image, int lo, int hi) { int j, i; if (logptr) { @@ -225,17 +217,18 @@ print_context(char *title, struct driz_param_t *p, int lo, int hi) { void print_pixmap(char *title, struct driz_param_t *p, int lo, int hi) { - int i, j, k; + int i, j, k; char *axis[2] = {"x", "y"}; if (logptr) { - for (k = 0; k < 2; k ++) { + for (k = 0; k < 2; k++) { fprintf(logptr, "\n%s %s axis\n", title, axis[k]); - for (j = 0; j < image_size[1]; j++ ) { + for (j = 0; j < image_size[1]; j++) { for (i = 0; i < image_size[0]; i++) { if (i >= lo && i < hi && j >= lo && j < hi) { - fprintf(logptr, "%10.2f", get_pixmap(p->pixmap, i, j)[k]); + fprintf(logptr, "%10.2f", + get_pixmap(p->pixmap, i, j)[k]); } } @@ -254,7 +247,7 @@ setup_parameters() { /* Initialize the parameter struct with vanilla defaults */ struct driz_param_t *p; - p = (struct driz_param_t *) malloc(sizeof(struct driz_param_t)); + p = (struct driz_param_t *)malloc(sizeof(struct driz_param_t)); driz_param_init(p); @@ -280,7 +273,7 @@ setup_parameters() { p->nmiss = 0; p->nskip = 0; - error = (struct driz_error_t *) malloc(sizeof(struct driz_error_t)); + error = (struct driz_error_t *)malloc(sizeof(struct driz_error_t)); driz_error_init(error); p->error = error; init_pixmap(p); @@ -302,7 +295,6 @@ setup_parameters() { void teardown_parameters(struct driz_param_t *p) { - if (logptr) { fclose(logptr); logptr = NULL; @@ -312,24 +304,17 @@ teardown_parameters(struct driz_param_t *p) { free(p); } -FCT_BGN_FN(utest_cdrizzle) -{ - FCT_FIXTURE_SUITE_BGN("unit tests for drizzle") - { - FCT_SETUP_BGN() - { - } +FCT_BGN_FN(utest_cdrizzle) { + FCT_FIXTURE_SUITE_BGN("unit tests for drizzle") { + FCT_SETUP_BGN() {} FCT_SETUP_END(); - FCT_TEARDOWN_BGN() - { - } + FCT_TEARDOWN_BGN() {} FCT_TEARDOWN_END(); - FCT_TEST_BGN(utest_shrink_bbox) - { + FCT_TEST_BGN(utest_shrink_bbox) { int xmin, xmax, ymin, ymax; - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); nan_pixmap(p); @@ -359,8 +344,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_lookup_01) - { + FCT_TEST_BGN(utest_map_lookup_01) { struct driz_param_t *p; p = setup_parameters(); stretch_pixmap(p, 1000.0); @@ -375,8 +359,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_lookup_02) - { + FCT_TEST_BGN(utest_map_lookup_02) { struct driz_param_t *p; p = setup_parameters(); stretch_pixmap(p, 1000.0); @@ -391,8 +374,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_lookup_03) - { + FCT_TEST_BGN(utest_map_lookup_03) { struct driz_param_t *p; p = setup_parameters(); stretch_pixmap(p, 1000.0); @@ -406,8 +388,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_lookup_04) - { + FCT_TEST_BGN(utest_map_lookup_04) { struct driz_param_t *p; p = setup_parameters(); stretch_pixmap(p, 1000.0); @@ -421,8 +402,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_point_01) - { + FCT_TEST_BGN(utest_map_point_01) { double ix, iy, ox, oy; struct driz_param_t *p; @@ -441,8 +421,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_point_02) - { + FCT_TEST_BGN(utest_map_point_02) { double ix, iy, ox, oy; struct driz_param_t *p; @@ -461,8 +440,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_point_03) - { + FCT_TEST_BGN(utest_map_point_03) { double ix, iy, ox, oy; int status; struct driz_param_t *p; @@ -490,8 +468,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_map_point_04) - { + FCT_TEST_BGN(utest_map_point_04) { double ix, iy, ox, oy; int status; struct driz_param_t *p; @@ -520,20 +497,19 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_check_line_overlap_01) - { + FCT_TEST_BGN(utest_check_line_overlap_01) { struct scanner s; int ymin, ymax, shift, status; /* Test for complete overlap */ - const integer_t j = 0; /* which image line to check ? */ - integer_t xbounds[2]; /* start of in-bounds */ - struct driz_param_t *p; /* parameter structure */ + const integer_t j = 0; /* which image line to check ? */ + integer_t xbounds[2]; /* start of in-bounds */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); shift = 0; - offset_pixmap(p, (double) shift, 0.0); + offset_pixmap(p, (double)shift, 0.0); init_image_scanner(p, &s, &ymin, &ymax); status = get_scanline_limits(&s, j, &xbounds[0], &xbounds[1]); @@ -545,25 +521,23 @@ FCT_BGN_FN(utest_cdrizzle) fct_chk_eq_int(xbounds[1], MAX(0, MIN(image_size[0] - 1 - shift, image_size[0] - 1))); - teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_check_line_overlap_02) - { + FCT_TEST_BGN(utest_check_line_overlap_02) { struct scanner s; int ymin, ymax, shift, status; /* Test for half overlap */ - const integer_t j = 0; /* which image line to check ? */ - integer_t xbounds[2]; /* start of in-bounds */ - struct driz_param_t *p; /* parameter structure */ + const integer_t j = 0; /* which image line to check ? */ + integer_t xbounds[2]; /* start of in-bounds */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); shift = 70; - offset_pixmap(p, (double) shift, 0.0); + offset_pixmap(p, (double)shift, 0.0); init_image_scanner(p, &s, &ymin, &ymax); status = get_scanline_limits(&s, j, &xbounds[0], &xbounds[1]); @@ -579,20 +553,19 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_check_line_overlap_03) - { + FCT_TEST_BGN(utest_check_line_overlap_03) { struct scanner s; int ymin, ymax, shift, status; /* Test for negative half overlap */ - const integer_t j = 0; /* which image line to check ? */ - integer_t xbounds[2]; /* start of in-bounds */ - struct driz_param_t *p; /* parameter structure */ + const integer_t j = 0; /* which image line to check ? */ + integer_t xbounds[2]; /* start of in-bounds */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); shift = -70; - offset_pixmap(p, (double) shift, 0.0); + offset_pixmap(p, (double)shift, 0.0); init_image_scanner(p, &s, &ymin, &ymax); status = get_scanline_limits(&s, j, &xbounds[0], &xbounds[1]); @@ -608,18 +581,17 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_check_image_overlap_02) - { + FCT_TEST_BGN(utest_check_image_overlap_02) { struct scanner s; int ymin, ymax, shift; /* Test for half overlap */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); shift = 70; - offset_pixmap(p, 0.0, (double) shift); + offset_pixmap(p, 0.0, (double)shift); init_image_scanner(p, &s, &ymin, &ymax); @@ -631,18 +603,17 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_check_image_overlap_03) - { + FCT_TEST_BGN(utest_check_image_overlap_03) { struct scanner s; int ymin, ymax, shift; /* Test for negative half overlap */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ p = setup_parameters(); shift = -70; - offset_pixmap(p, 0.0, (double) shift); + offset_pixmap(p, 0.0, (double)shift); init_image_scanner(p, &s, &ymin, &ymax); @@ -654,8 +625,7 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_compute_area_01) - { + FCT_TEST_BGN(utest_compute_area_01) { /* Test compute area with aligned square entirely inside */ double area; double is, js, x[4], y[4]; @@ -674,12 +644,10 @@ FCT_BGN_FN(utest_cdrizzle) area = compute_area(is, js, x, y); fct_chk_eq_dbl(area, 0.25); - } FCT_TEST_END(); - FCT_TEST_BGN(utest_compute_area_02) - { + FCT_TEST_BGN(utest_compute_area_02) { /* Test compute area with diagonal square entirely inside */ double area; double is, js, x[4], y[4]; @@ -698,12 +666,10 @@ FCT_BGN_FN(utest_cdrizzle) area = compute_area(is, js, x, y); fct_chk_eq_dbl(area, 0.125); - } FCT_TEST_END(); - FCT_TEST_BGN(utest_compute_area_03) - { + FCT_TEST_BGN(utest_compute_area_03) { /* Test compute area with aligned square with overlap */ double area; double is, js, x[4], y[4]; @@ -722,12 +688,10 @@ FCT_BGN_FN(utest_cdrizzle) area = compute_area(is, js, x, y); fct_chk_eq_dbl(area, 0.25); - } FCT_TEST_END(); - FCT_TEST_BGN(utest_compute_area_04) - { + FCT_TEST_BGN(utest_compute_area_04) { /* Test compute area with diagonal square with overlap */ double area; double is, js, x[4], y[4]; @@ -746,37 +710,41 @@ FCT_BGN_FN(utest_cdrizzle) area = compute_area(is, js, x, y); fct_chk_eq_dbl(area, 0.0625); - } FCT_TEST_END(); - FCT_TEST_BGN(utest_compute_area_05) - { + FCT_TEST_BGN(utest_compute_area_05) { /* Test compute area with marching diagonal square */ int i, j; double is, js, x[4], y[4], area; - double area_ok[7][7] = - {{0.125000, 0.218750, 0.250000, 0.218750, 0.125000, 0.031250, 0.000000}, - {0.218750, 0.375000, 0.437500, 0.375000, 0.218750, 0.062500, 0.000000}, - {0.250000, 0.437500, 0.500000, 0.437500, 0.250000, 0.062500, 0.000000}, - {0.218750, 0.375000, 0.437500, 0.375000, 0.218750, 0.062500, 0.000000}, - {0.125000, 0.218750, 0.250000, 0.218750, 0.125000, 0.031250, 0.000000}, - {0.031250, 0.062500, 0.062500, 0.062500, 0.031250, 0.000000, 0.000000}, - {0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000}}; + double area_ok[7][7] = {{0.125000, 0.218750, 0.250000, 0.218750, + 0.125000, 0.031250, 0.000000}, + {0.218750, 0.375000, 0.437500, 0.375000, + 0.218750, 0.062500, 0.000000}, + {0.250000, 0.437500, 0.500000, 0.437500, + 0.250000, 0.062500, 0.000000}, + {0.218750, 0.375000, 0.437500, 0.375000, + 0.218750, 0.062500, 0.000000}, + {0.125000, 0.218750, 0.250000, 0.218750, + 0.125000, 0.031250, 0.000000}, + {0.031250, 0.062500, 0.062500, 0.062500, + 0.031250, 0.000000, 0.000000}, + {0.000000, 0.000000, 0.000000, 0.000000, + 0.000000, 0.000000, 0.000000}}; is = 1.0; js = 1.0; - for (i = 0; i <= 6; ++ i) { - for (j = 0; j <= 6; ++ j) { - x[0] = 0.25 * (double) i; - y[0] = 0.25 * (double) j + 0.5; - x[1] = 0.25 * (double) i + 0.5; - y[1] = 0.25 * (double) j; - x[2] = 0.25 * (double) i + 1.0; - y[2] = 0.25 * (double) j + 0.5; - x[3] = 0.25 * (double) i + 0.5; - y[3] = 0.25 * (double) j + 1.0; + for (i = 0; i <= 6; ++i) { + for (j = 0; j <= 6; ++j) { + x[0] = 0.25 * (double)i; + y[0] = 0.25 * (double)j + 0.5; + x[1] = 0.25 * (double)i + 0.5; + y[1] = 0.25 * (double)j; + x[2] = 0.25 * (double)i + 1.0; + y[2] = 0.25 * (double)j + 0.5; + x[3] = 0.25 * (double)i + 0.5; + y[3] = 0.25 * (double)j + 1.0; area = compute_area(is, js, x, y); @@ -787,13 +755,12 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_01) - { + FCT_TEST_BGN(utest_do_kernel_square_01) { /* Simplest case */ - integer_t x1; /* start of in-bounds */ - integer_t j, x2; /* end of in-bounds */ - struct driz_param_t *p; /* parameter structure */ + integer_t x1; /* start of in-bounds */ + integer_t j, x2; /* end of in-bounds */ + struct driz_param_t *p; /* parameter structure */ int n, status; n = 100; @@ -808,18 +775,19 @@ FCT_BGN_FN(utest_cdrizzle) x2 = n; fct_chk_eq_int(status, 0); - fct_chk_eq_dbl(get_pixel(p->output_data, x1, j), get_pixel(p->data, x1, j)); - fct_chk_eq_dbl(get_pixel(p->output_data, x2-1, j), get_pixel(p->data, x2-1, j)); + fct_chk_eq_dbl(get_pixel(p->output_data, x1, j), + get_pixel(p->data, x1, j)); + fct_chk_eq_dbl(get_pixel(p->output_data, x2 - 1, j), + get_pixel(p->data, x2 - 1, j)); teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_02) - { + FCT_TEST_BGN(utest_do_kernel_square_02) { /* Offset image */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ integer_t j; int k, n, status; double offset; @@ -834,19 +802,19 @@ FCT_BGN_FN(utest_cdrizzle) status = do_kernel_square(p); fct_chk_eq_int(status, 0); - for (k = 1; k < n-2; k++) { - fct_chk_eq_dbl(get_pixel(p->output_data, (k+1), j), get_pixel(p->data, k, j)); + for (k = 1; k < n - 2; k++) { + fct_chk_eq_dbl(get_pixel(p->output_data, (k + 1), j), + get_pixel(p->data, k, j)); } teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_03) - { + FCT_TEST_BGN(utest_do_kernel_square_03) { /* Single pixel set */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ integer_t j; int k, n; double offset, value; @@ -864,17 +832,17 @@ FCT_BGN_FN(utest_cdrizzle) k = 3; do_kernel_square(p); - fct_chk_eq_dbl(get_pixel(p->output_data, (k+2), k), get_pixel(p->data, k, k)); + fct_chk_eq_dbl(get_pixel(p->output_data, (k + 2), k), + get_pixel(p->data, k, k)); teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_04) - { + FCT_TEST_BGN(utest_do_kernel_square_04) { /* Single pixel, fractional pixel offset */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ integer_t i, j, k; double offset, value; @@ -890,7 +858,8 @@ FCT_BGN_FN(utest_cdrizzle) for (i = 2; i <= 3; ++i) { for (j = 2; j <= 3; ++j) { - fct_chk_eq_dbl(get_pixel(p->output_data, (i+k), (j+k)), value/4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, (i + k), (j + k)), + value / 4.0); } } @@ -898,12 +867,11 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_05) - { + FCT_TEST_BGN(utest_do_kernel_square_05) { /* Diagonal line, fractional pixel offset */ - struct driz_param_t *p; /* parameter structure */ - integer_t i,j; + struct driz_param_t *p; /* parameter structure */ + integer_t i, j; int n; double offset, value; @@ -913,28 +881,29 @@ FCT_BGN_FN(utest_cdrizzle) p = setup_parameters(); offset_pixmap(p, offset, offset); - for (j = 1; j < n-1; ++j) { + for (j = 1; j < n - 1; ++j) { set_pixel(p->data, j, j, value); } do_kernel_square(p); for (i = 4; i < n - 2; ++i) { - fct_chk_eq_dbl(get_pixel(p->output_data, i, i), value/2.0); - fct_chk_eq_dbl(get_pixel(p->output_data, i-1, i), value/4.0); - fct_chk_eq_dbl(get_pixel(p->output_data, i, i-1), value/4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i, i), value / 2.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i - 1, i), + value / 4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i, i - 1), + value / 4.0); } teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_06) - { + FCT_TEST_BGN(utest_do_kernel_square_06) { /* Block of pixels, whole number offset */ - struct driz_param_t *p; /* parameter structure */ - integer_t i,j; + struct driz_param_t *p; /* parameter structure */ + integer_t i, j; int k; double offset, value; @@ -954,7 +923,8 @@ FCT_BGN_FN(utest_cdrizzle) for (i = 0; i < k; ++i) { for (j = 0; j < k; ++j) { - fct_chk_eq_dbl(get_pixel(p->output_data, (i+2), (j+2)), value); + fct_chk_eq_dbl(get_pixel(p->output_data, (i + 2), (j + 2)), + value); } } @@ -962,11 +932,10 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_do_kernel_square_07) - { + FCT_TEST_BGN(utest_do_kernel_square_07) { /* Block of pixels, fractional offset */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ integer_t i, j; int k; double offset, value; @@ -977,8 +946,8 @@ FCT_BGN_FN(utest_cdrizzle) p = setup_parameters(); offset_pixmap(p, offset, offset); - for (i = 1; i < k+1; ++i) { - for (j = 1; j < k+1; ++j) { + for (i = 1; i < k + 1; ++i) { + for (j = 1; j < k + 1; ++j) { set_pixel(p->data, i, j, value); } } @@ -999,11 +968,10 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_dobox_01) - { + FCT_TEST_BGN(utest_dobox_01) { /* Single pixel set, whole number offset */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int k; double offset, value; @@ -1018,16 +986,16 @@ FCT_BGN_FN(utest_cdrizzle) set_pixel(p->data, k, k, value); dobox(p); - fct_chk_eq_dbl(get_pixel(p->output_data, (k+2), (k+2)), get_pixel(p->data, k, k)); + fct_chk_eq_dbl(get_pixel(p->output_data, (k + 2), (k + 2)), + get_pixel(p->data, k, k)); fct_chk_eq_int(p->nskip, 2); } FCT_TEST_END(); - FCT_TEST_BGN(utest_dobox_02) - { + FCT_TEST_BGN(utest_dobox_02) { /* Single pixel, fractional pixel offset */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int i, j, k; double offset, value; @@ -1044,17 +1012,17 @@ FCT_BGN_FN(utest_cdrizzle) for (i = 2; i <= 3; ++i) { for (j = 2; j <= 3; ++j) { - fct_chk_eq_dbl(get_pixel(p->output_data, (i+k), (j+k)), value/4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, (i + k), (j + k)), + value / 4.0); } } } FCT_TEST_END(); - FCT_TEST_BGN(utest_dobox_03) - { + FCT_TEST_BGN(utest_dobox_03) { /* Turbo mode kernel, diagonal line of pixels set */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int i, j, n; double offset, value; @@ -1073,20 +1041,21 @@ FCT_BGN_FN(utest_cdrizzle) dobox(p); for (i = 4; i < n; ++i) { - fct_chk_eq_dbl(get_pixel(p->output_data, i, i), value/2.0); - fct_chk_eq_dbl(get_pixel(p->output_data, i-1, i), value/4.0); - fct_chk_eq_dbl(get_pixel(p->output_data, i, i-1), value/4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i, i), value / 2.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i - 1, i), + value / 4.0); + fct_chk_eq_dbl(get_pixel(p->output_data, i, i - 1), + value / 4.0); } teardown_parameters(p); } FCT_TEST_END(); - FCT_TEST_BGN(utest_dobox_04) - { + FCT_TEST_BGN(utest_dobox_04) { /* Check that context map is set for the affected pixels */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int i, j, n; double offset, value; integer_t bv; @@ -1099,7 +1068,7 @@ FCT_BGN_FN(utest_cdrizzle) offset_pixmap(p, offset, offset); // DBG p->kernel = kernel_turbo; - for (j = 1; j < n-1; ++j) { + for (j = 1; j < n - 1; ++j) { set_pixel(p->data, j, j, value); } @@ -1113,11 +1082,10 @@ FCT_BGN_FN(utest_cdrizzle) } FCT_TEST_END(); - FCT_TEST_BGN(utest_doblot_01) - { + FCT_TEST_BGN(utest_doblot_01) { /* Single pixel set blinear interpolation */ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int k; double offset, value; @@ -1133,15 +1101,15 @@ FCT_BGN_FN(utest_cdrizzle) set_pixel(p->data, k, k, value); doblot(p); - fct_chk_eq_dbl(get_pixel(p->output_data, (k+2), (k+2)), get_pixel(p->data, k, k)); + fct_chk_eq_dbl(get_pixel(p->output_data, (k + 2), (k + 2)), + get_pixel(p->data, k, k)); } FCT_TEST_END(); - FCT_TEST_BGN(utest_doblot_02) - { + FCT_TEST_BGN(utest_doblot_02) { /* Single pixel set quintic interpolation*/ - struct driz_param_t *p; /* parameter structure */ + struct driz_param_t *p; /* parameter structure */ int k; double offset, value; @@ -1157,11 +1125,11 @@ FCT_BGN_FN(utest_cdrizzle) set_pixel(p->data, k, k, value); doblot(p); - fct_chk_eq_dbl(get_pixel(p->output_data, (k+2), (k+2)), get_pixel(p->data, k, k)); + fct_chk_eq_dbl(get_pixel(p->output_data, (k + 2), (k + 2)), + get_pixel(p->data, k, k)); } FCT_TEST_END(); - - } + } FCT_FIXTURE_SUITE_END(); } FCT_END_FN();