/* $Id: cpl_tools.c,v 1.96 2008/01/28 09:33:42 yjung Exp $ * * This file is part of the ESO Common Pipeline Library * Copyright (C) 2001-2004 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: yjung $ * $Date: 2008/01/28 09:33:42 $ * $Revision: 1.96 $ * $Name: $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes -----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include #include #include "time.h" /* Needed by strcmp(), strncmp(), strlen() */ #include #include #include #include #include "cpl_propertylist_impl.h" #include "cpl_memory.h" #include "cpl_tools.h" /*----------------------------------------------------------------------------- Defines -----------------------------------------------------------------------------*/ /* Swap macros */ #define CPL_DOUBLE_SWAP(a,b) { register double t=(a);(a)=(b);(b)=t; } #define CPL_FLOAT_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; } #define CPL_INT_SWAP(a,b) { register int t=(a);(a)=(b);(b)=t; } /*----------------------------------------------------------------------------*/ /** @enum cpl_fits_keytype @internal @brief Possible FITS key types This enum stores all possible types for a FITS keyword. These determine the order of appearance in a header, they are a crucial point for DICB (ESO) compatibility. This classification is internal to this module. */ /*----------------------------------------------------------------------------*/ typedef enum _cpl_fits_keytype_ { cpl_fits_keytype_undef =0, cpl_fits_keytype_top =1, /* Mandatory keywords */ /* All FITS files */ cpl_fits_keytype_bitpix =2, cpl_fits_keytype_naxis =3, cpl_fits_keytype_naxis1 =11, cpl_fits_keytype_naxis2 =12, cpl_fits_keytype_naxis3 =13, cpl_fits_keytype_naxis4 =14, cpl_fits_keytype_naxisi =20, /* Random groups only */ cpl_fits_keytype_group =30, /* Extensions */ cpl_fits_keytype_pcount =31, cpl_fits_keytype_gcount =32, /* Main header */ cpl_fits_keytype_extend =33, /* Images */ cpl_fits_keytype_bscale =34, cpl_fits_keytype_bzero =35, /* Tables */ cpl_fits_keytype_tfields =36, cpl_fits_keytype_tbcoli =40, cpl_fits_keytype_tformi =41, /* Other primary keywords */ cpl_fits_keytype_primary =100, /* HIERARCH ESO keywords ordered according to DICB */ cpl_fits_keytype_hierarch_dpr =200, cpl_fits_keytype_hierarch_obs =201, cpl_fits_keytype_hierarch_tpl =202, cpl_fits_keytype_hierarch_gen =203, cpl_fits_keytype_hierarch_tel =204, cpl_fits_keytype_hierarch_ins =205, cpl_fits_keytype_hierarch_det =206, cpl_fits_keytype_hierarch_log =207, cpl_fits_keytype_hierarch_pro =208, /* Other HIERARCH keywords */ cpl_fits_keytype_hierarch =300, /* HISTORY and COMMENT */ cpl_fits_keytype_history =400, cpl_fits_keytype_comment =500, /* END */ cpl_fits_keytype_end =1000 } cpl_fits_keytype ; /*----------------------------------------------------------------------------- Private function prototypes -----------------------------------------------------------------------------*/ static double cpl_tools_get_median_9double(double *); static float cpl_tools_get_median_9float(float *); static int cpl_fits_property_comparison(const cpl_property *, const cpl_property *) ; static cpl_fits_keytype cpl_fits_property_get_type(const char *) ; /*----------------------------------------------------------------------------- Private variables -----------------------------------------------------------------------------*/ static cpl_flops flop_count = 0; /*----------------------------------------------------------------------------*/ /** * @defgroup cpl_tools Utility functions * * This module provides various functions to apply simple operations on * data arrays (sorting, median computation) and a function to monitor the * CPU time. * * @par Synopsis: * @code * #include "cpl_tools.h" * @endcode */ /*----------------------------------------------------------------------------*/ /**@{*/ /*----------------------------------------------------------------------------- Function codes -----------------------------------------------------------------------------*/ /*----------------------------------------------------------------------------*/ /** @internal @brief Get the kth smallest value in a double array @param a the double array @param n the array size @param k the requested value position in the sorted array @return the kth smallest value in the sorted array. The input array is modified. In case of error, the #_cpl_error_code_ code is set, and the returned value is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT */ /*----------------------------------------------------------------------------*/ double cpl_tools_get_kth_double( double * a, int n, int k) { register double x ; register int i, j, l, m ; cpl_ensure(a, CPL_ERROR_NULL_INPUT, 0.00) ; l=0 ; m=n-1 ; while (l=1 ; i--) { if (pix_arr[i-1] <= a) break; pix_arr[i] = pix_arr[i-1]; } pix_arr[i] = a; } if (j_stack == 0) break; ir = i_stack[j_stack-- -1]; l = i_stack[j_stack-- -1]; } else { k = (l+ir) >> 1; CPL_DOUBLE_SWAP(pix_arr[k-1], pix_arr[l]) if (pix_arr[l] > pix_arr[ir-1]) { CPL_DOUBLE_SWAP(pix_arr[l], pix_arr[ir-1]) } if (pix_arr[l-1] > pix_arr[ir-1]) { CPL_DOUBLE_SWAP(pix_arr[l-1], pix_arr[ir-1]) } if (pix_arr[l] > pix_arr[l-1]) { CPL_DOUBLE_SWAP(pix_arr[l], pix_arr[l-1]) } i = l+1; j = ir; a = pix_arr[l-1]; for (;;) { do i++; while (pix_arr[i-1] < a); do j--; while (pix_arr[j-1] > a); if (j < i) break; CPL_DOUBLE_SWAP(pix_arr[i-1], pix_arr[j-1]); } pix_arr[l-1] = pix_arr[j-1]; pix_arr[j-1] = a; j_stack += 2; if (j_stack > CPL_PIX_STACK_SIZE) { /* Should never reach here */ return CPL_ERROR_ILLEGAL_INPUT ; } if (ir-i+1 >= j-l) { i_stack[j_stack-1] = ir; i_stack[j_stack-2] = i; ir = j-1; } else { i_stack[j_stack-1] = j-1; i_stack[j_stack-2] = l; l = i; } } } return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @internal @brief Sort an integer array @param pix_arr the array to sort @param n the array size @return the #_cpl_error_code_ or CPL_ERROR_NONE Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT - CPL_ERROR_ILLEGAL_INPUT */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_tools_sort_int( int * pix_arr, int n) { int i, ir, j, k, l; int i_stack[CPL_PIX_STACK_SIZE] ; int j_stack ; int a ; /* Check entries */ cpl_ensure(pix_arr, CPL_ERROR_NULL_INPUT, CPL_ERROR_NULL_INPUT) ; ir = n ; l = 1 ; j_stack = 0 ; for (;;) { if (ir-l < 7) { for (j=l+1 ; j<=ir ; j++) { a = pix_arr[j-1]; for (i=j-1 ; i>=1 ; i--) { if (pix_arr[i-1] <= a) break; pix_arr[i] = pix_arr[i-1]; } pix_arr[i] = a; } if (j_stack == 0) break; ir = i_stack[j_stack-- -1]; l = i_stack[j_stack-- -1]; } else { k = (l+ir) >> 1; CPL_INT_SWAP(pix_arr[k-1], pix_arr[l]) if (pix_arr[l] > pix_arr[ir-1]) { CPL_INT_SWAP(pix_arr[l], pix_arr[ir-1]) } if (pix_arr[l-1] > pix_arr[ir-1]) { CPL_INT_SWAP(pix_arr[l-1], pix_arr[ir-1]) } if (pix_arr[l] > pix_arr[l-1]) { CPL_INT_SWAP(pix_arr[l], pix_arr[l-1]) } i = l+1; j = ir; a = pix_arr[l-1]; for (;;) { do i++; while (pix_arr[i-1] < a); do j--; while (pix_arr[j-1] > a); if (j < i) break; CPL_INT_SWAP(pix_arr[i-1], pix_arr[j-1]); } pix_arr[l-1] = pix_arr[j-1]; pix_arr[j-1] = a; j_stack += 2; if (j_stack > CPL_PIX_STACK_SIZE) { /* Should never reach here */ return CPL_ERROR_ILLEGAL_INPUT ; } if (ir-i+1 >= j-l) { i_stack[j_stack-1] = ir; i_stack[j_stack-2] = i; ir = j-1; } else { i_stack[j_stack-1] = j-1; i_stack[j_stack-2] = l; l = i; } } } return CPL_ERROR_NONE ; } #undef CPL_PIX_STACK_SIZE /*----------------------------------------------------------------------------*/ /** @internal @brief Compute the sample variance of a double array @param a The double array @param n The non-negative array size @param mean The sample mean of a @return The variance The sample variance is S(n) = (1/n) sum((x_i-mean)^2) (i=1 -> n) mean may be different from the actual mean, in which case something else is computed. */ /*----------------------------------------------------------------------------*/ double cpl_tools_get_variance_double(const double * a, int n, double mean) { double variance = 0; int i; cpl_ensure(a != NULL, CPL_ERROR_NULL_INPUT, 0.0); cpl_ensure(n >= 0, CPL_ERROR_ILLEGAL_INPUT, 0.0); /* Compute the variance of the elements using the recurrence relation var_n = var_{n-1} + ((x_n-mean)*(x_n-mean) - var_{n-1})/(n+1) */ for (i=0; i < n; i++) { const double delta = a[i] - mean; variance += (delta * delta - variance) / (i + 1); } return variance; } /*----------------------------------------------------------------------------*/ /** @internal @brief Get the median value in a double array @param a the double array @param n the array size @return the requested median @see cpl_tools_get_kth_double() For a finite population or sample, the median is the middle value of an odd number of values (arranged in ascending order) or any value between the two middle values of an even number of values. The computed median should actually be a value of the input array, so, in case of an even number of values in the input array, the median would be the lowest of the two central values. The function is optimized for 9 elements. Example: the median of (1 2 3) is 2 and the median of (1 2 3 4) is 2. */ /*----------------------------------------------------------------------------*/ double cpl_tools_get_median_double( double * a, int n) { return n == 9 ? cpl_tools_get_median_9double(a) : cpl_tools_get_kth_double(a, n, ((n&1) ? (n/2) : ((n/2)-1))); } /*----------------------------------------------------------------------------*/ /** @internal @brief Get the median value in a float array @param a the float array @param n the array size @return the requested median @see cpl_tools_get_median_double() */ /*----------------------------------------------------------------------------*/ float cpl_tools_get_median_float( float * a, int n) { return n == 9 ? cpl_tools_get_median_9float(a) : cpl_tools_get_kth_float(a, n, ((n&1) ? (n/2) : ((n/2)-1))); } /*----------------------------------------------------------------------------*/ /** @internal @brief Get the median value in an integer array @param a the integer array @param n the array size @return the requested median @see cpl_tools_get_median_double() */ /*----------------------------------------------------------------------------*/ int cpl_tools_get_median_int( int * a, int n) { return cpl_tools_get_kth_int(a,n,(((n)&1)?((n)/2):(((n)/2)-1))) ; } /*----------------------------------------------------------------------------*/ /** @internal @brief Timer handling, including benchmarking. @param mode [CPL_CLOCK_START CPL_CLOCK_STOP CPL_CLOCK_TIME CPL_CLOCK_RESET] @return An elapsed time in seconds as a double. This routine is useful to benchmark the time spent in an image processing routine. It is used in three steps: 1) Start the timer by calling @c cpl_tools_get_cputime(CPL_CLOCK_START) 2) Stop the timer by calling @c cpl_tools_get_cputime(CPL_CLOCK_STOP) 3) Get the total time by calling @c cpl_tools_get_cputime(CPL_CLOCK_TIME) Steps 1) and 2) may be repeated pairwise more than once. After any instance of step 2) step 3) may be repeated more than once. cpl_tools_get_cputime(CPL_CLOCK_TIME) may also be called before step 1), in which case it returns 0. The call with CPL_CLOCK_START returns the total elapsed time in seconds since the begining. The call with CPL_CLOCK_STOP returns the time elapsed since the last CPL_CLOCK_START call. The call with CPL_CLOCK_TIME return the total time elapsed between all previous CPL_CLOCK_START - CPL_CLOCK_STOP calls. Additionally, the call with CPL_CLOCK_RESET sets the internal state of the timing function to its initial state, something that may be necessary when recovering from an error. In this case the function returns zero. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_ILLEGAL_INPUT */ /*----------------------------------------------------------------------------*/ double cpl_tools_get_cputime(cpl_clock mode) { /* The state of the static variables */ static cpl_clock pmode = CPL_CLOCK_STOP; /* clock() may wrap around after as little as 36 minutes, so store the result in a double */ static clock_t chrono = 0; static double total = 0; clock_t diff; double result ; switch (mode) { case CPL_CLOCK_START: cpl_ensure(pmode != CPL_CLOCK_START, CPL_ERROR_ILLEGAL_INPUT, 0); chrono = clock() ; result = chrono; pmode = mode; break ; case CPL_CLOCK_STOP: cpl_ensure(pmode == CPL_CLOCK_START, CPL_ERROR_ILLEGAL_INPUT, 0); /* At this point clock() may have wrapped around, causing it to be much smaller than chrono. Subtraction must therefore be done in type (clock_t). */ diff = clock() - chrono; result = diff; total += result; pmode = mode; break ; case CPL_CLOCK_TIME: cpl_ensure(pmode != CPL_CLOCK_START, CPL_ERROR_ILLEGAL_INPUT, 0); result = total ; break ; case CPL_CLOCK_RESET: /* Set the static variables to their initial values */ pmode = CPL_CLOCK_STOP; chrono = 0; total = 0.0; result = 0.0; /* Need to return something */ break ; default: /* Only reached if cpl_clock is extended in error */ cpl_ensure(0, CPL_ERROR_UNSUPPORTED_MODE, 0.0); } return result / CLOCKS_PER_SEC; } /*----------------------------------------------------------------------------*/ /** @internal @brief Determine if an integer is a power of 2. @param p Integer to check. @return The corresponding power of 2 or -1 if p is not a power of 2 */ /*----------------------------------------------------------------------------*/ int cpl_tools_is_power_of_2(int p) { int ipow = 0; int done; if (p <= 0) return -1; /* Count until 1st non-zero bit is found */ do { done = p & 1; p >>= 1; ipow++; } while (!done); /* The number is a power iff p is now zero */ return p == 0 ? ipow-1 : -1; } /*----------------------------------------------------------------------------*/ /** @internal @brief Compute x to the power of y @param x The base @param y The power @return x to the power of y @note Iff y is negative and x is zero an actual division by zero will occur Apart from a possible difference in round-off the result equals pow(x,y). */ /*----------------------------------------------------------------------------*/ double cpl_tools_ipow(double x, int y) { double result; double pow2 = x; int p = abs(y); /* Compute the result as a product of powers of 2 of x. - this method may produce (slightly) different results than pow(x,y) */ /* Handle least significant bit in p here in order to avoid an unnecessary multiplication of pow2 - which could cause an over- or underflow */ /* Also, 0^0 is 1, while any other power of 0 is 0 */ result = p & 1 ? x : 1.0; while (p >>= 1) { pow2 *= pow2; /* Process least significant bit in p */ if (p & 1) result *= pow2; } /* It is left to the caller to ensure that no division by zero occurs */ return y < 0 ? 1.0/result : result; } /*----------------------------------------------------------------------------*/ /** @internal @brief Return a string with a CFITSIO error message @param status A CFITSIO status variable @return A static string with the message @note It is forbidden to modify or free this string */ /*----------------------------------------------------------------------------*/ const char * cpl_tools_get_cfitsio_msg(int status) { static char cfitsio_msg[FLEN_STATUS]; fits_get_errstatus(status, cfitsio_msg); /* The CFITSIO documentation does not guarantee that the string is properly terminated */ /* Actually, a look at the fits_get_errstatus() source code (versions 2.510 and 3.030) reveals that the error message is written to the caller- supplied buffer with strcpy(). So the string will be NULL-terminated - and the CFITSIO developers better remember to keep FLEN_STATUS big enough for their messages, or they will cause a segmentation violation. */ cfitsio_msg[FLEN_STATUS-1] = '\0'; /* Better safe than sorry */ return cfitsio_msg; } /*----------------------------------------------------------------------------*/ /** @internal @brief Generate a valid FITS header for saving functions @param fptr The fitsfile to modify @param keywords The set of keywords to add to the minimal header @param to_rm The set of keys to remove @return 1 newly allocated valid header The passed file should contain a minimal header. The propertylist is sorted (DICB) and added after the miniaml header. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if fptr is NULL - CPL_ERROR_ILLEGAL_INPUT if the propertylist cannot be written */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_fits_add_properties( fitsfile * fptr, const cpl_propertylist * keywords, const char * to_rm) { cpl_propertylist * out ; /* Test entries */ if (keywords == NULL) return CPL_ERROR_NONE ; cpl_ensure_code(fptr, CPL_ERROR_NULL_INPUT) ; /* * The following has been replaced by the use of the new * cpl_propertylist_to_fitsfile(), which takes a to_rm * as an argument. */ /* Copy all but the black-listed properties */ /* if (to_rm != NULL) { out = cpl_propertylist_new(); if (cpl_propertylist_copy_property_regexp(out, keywords, to_rm, 1)) { cpl_propertylist_delete(out); cpl_ensure_code(0, cpl_error_get_code()); } } else { out = cpl_propertylist_duplicate(keywords) ; }*/ out = cpl_propertylist_duplicate(keywords); /* Sort the propertylist */ cpl_propertylist_sort(out, cpl_fits_property_comparison) ; /* Write the propertylist in the file */ if (cpl_propertylist_to_fitsfile(fptr, out, to_rm) != CPL_ERROR_NONE) { cpl_propertylist_delete(out) ; cpl_ensure_code(0, CPL_ERROR_ILLEGAL_INPUT) ; } cpl_propertylist_delete(out) ; return CPL_ERROR_NONE ; } /*----------------------------------------------------------------------------*/ /** @internal @brief Comparison function of two properties @param p1 the first property @param p2 the second property @return -1 if p1p2 This function implements the sorting of the FITS header keywords to insure DICB compliant products */ /*----------------------------------------------------------------------------*/ static int cpl_fits_property_comparison( const cpl_property * p1, const cpl_property * p2) { const char * key1 ; const char * key2 ; cpl_fits_keytype kt1, kt2 ; /* Get the keys */ key1 = cpl_property_get_name(p1) ; key2 = cpl_property_get_name(p2) ; /* Get the keywords types */ kt1 = cpl_fits_property_get_type(key1) ; kt2 = cpl_fits_property_get_type(key2) ; /* Compare the types */ if (kt1 > kt2) return 1 ; else if (kt1 < kt2) return -1 ; else return 0 ; } /*----------------------------------------------------------------------------*/ /** @internal @brief Get the FITS keyword type @param key the keyword @return the keyword type or -1 */ /*----------------------------------------------------------------------------*/ static cpl_fits_keytype cpl_fits_property_get_type(const char * key) { cpl_fits_keytype kt ; /* Check entries */ if (key == NULL) return cpl_fits_keytype_undef ; kt = cpl_fits_keytype_undef ; if (!strcmp(key, "SIMPLE") || !strcmp(key, "XTENSION")) kt = cpl_fits_keytype_top ; else if (!strcmp(key, "END")) kt = cpl_fits_keytype_end ; else if (!strcmp(key, "BITPIX")) kt = cpl_fits_keytype_bitpix ; else if (!strcmp(key, "NAXIS")) kt = cpl_fits_keytype_naxis ; else if (!strcmp(key, "NAXIS1")) kt = cpl_fits_keytype_naxis1 ; else if (!strcmp(key, "NAXIS2")) kt = cpl_fits_keytype_naxis2 ; else if (!strcmp(key, "NAXIS3")) kt = cpl_fits_keytype_naxis3 ; else if (!strcmp(key, "NAXIS4")) kt = cpl_fits_keytype_naxis4 ; else if (!strncmp(key, "NAXIS", 5)) kt = cpl_fits_keytype_naxisi ; else if (!strcmp(key, "GROUP")) kt = cpl_fits_keytype_group ; else if (!strcmp(key, "PCOUNT")) kt = cpl_fits_keytype_pcount ; else if (!strcmp(key, "GCOUNT")) kt = cpl_fits_keytype_gcount ; else if (!strcmp(key, "EXTEND")) kt = cpl_fits_keytype_extend ; else if (!strcmp(key, "BSCALE")) kt = cpl_fits_keytype_bscale ; else if (!strcmp(key, "BZERO")) kt = cpl_fits_keytype_bzero ; else if (!strcmp(key, "TFIELDS")) kt = cpl_fits_keytype_tfields ; else if (!strncmp(key, "TBCOL", 5)) kt = cpl_fits_keytype_tbcoli ; else if (!strncmp(key, "TFORM", 5)) kt = cpl_fits_keytype_tformi ; else if (!strncmp(key, "ESO DPR", 7)) kt = cpl_fits_keytype_hierarch_dpr ; else if (!strncmp(key, "ESO OBS", 7)) kt = cpl_fits_keytype_hierarch_obs ; else if (!strncmp(key, "ESO TPL", 7)) kt = cpl_fits_keytype_hierarch_tpl ; else if (!strncmp(key, "ESO GEN", 7)) kt = cpl_fits_keytype_hierarch_gen ; else if (!strncmp(key, "ESO TEL", 7)) kt = cpl_fits_keytype_hierarch_tel ; else if (!strncmp(key, "ESO INS", 7)) kt = cpl_fits_keytype_hierarch_ins ; else if (!strncmp(key, "ESO DET", 7)) kt = cpl_fits_keytype_hierarch_det; else if (!strncmp(key, "ESO LOG", 7)) kt = cpl_fits_keytype_hierarch_log ; else if (!strncmp(key, "ESO PRO", 7)) kt = cpl_fits_keytype_hierarch_pro ; else if (!strncmp(key, "ESO", 3)) kt = cpl_fits_keytype_hierarch ; else if (!strcmp(key, "HISTORY")) kt = cpl_fits_keytype_history ; else if (!strcmp(key, "COMMENT")) kt = cpl_fits_keytype_comment ; else if ((int)strlen(key)<9) kt = cpl_fits_keytype_primary ; return kt ; } /*----------------------------------------------------------------------------*/ /** @internal @brief Set the FLOPS counter to zero. @return void @note This function is intended to be used only by the CPL test modules. */ /*----------------------------------------------------------------------------*/ void cpl_tools_reset_flops (void) { flop_count = 0; } /*----------------------------------------------------------------------------*/ /** @internal @brief Add to the FLOPS counter. @param flops FLOPS to add @return void @note This function is intended to be used only by the CPL test modules. */ /*----------------------------------------------------------------------------*/ void cpl_tools_add_flops(cpl_flops flops) { flop_count += flops; } /*----------------------------------------------------------------------------*/ /** @internal @brief Get the FLOPS count. @return The FLOPS count @note This function is intended to be used only by the CPL test modules. */ /*----------------------------------------------------------------------------*/ cpl_flops cpl_tools_get_flops(void) { return flop_count; } /*----------------------------------------------------------------------------*/ /** @internal @brief Given p and u, modify the polynomial to p(x) := p(x+u) @param p The polynomial coefficients to be modified in place @param n The number of coefficients @param u The shift @return void @see cpl_polynomial_shift_1d @note The function will cx_assert() on NULL input. */ /*----------------------------------------------------------------------------*/ void cpl_polynomial_shift_double(double * coeffs, int n, double u) { int i, j; cx_assert( coeffs != NULL ); cx_assert( n > 0 ); for (j = 0; j < n-1; j++) for (i = 1; i < n - j; i++ ) coeffs[n-1-i] += coeffs[n-i] * u; cpl_tools_add_flops( n * ( n - 1) ); } /*----------------------------------------------------------------------------*/ /** @internal @brief Transform: xhat = x- mean(x) @param x The vector to be transformed @param pm On return, *pm is the mean of x @return The created transformed vector @note The function will cx_assert() on NULL input. */ /*----------------------------------------------------------------------------*/ cpl_vector * cpl_vector_transform_mean(const cpl_vector * x, double * pm) { cpl_vector * xhat = cpl_vector_duplicate(x); cx_assert( xhat != NULL ); cx_assert( pm != NULL ); *pm = cpl_vector_get_mean(xhat); cpl_vector_subtract_scalar(xhat, *pm); return xhat; } /*----------------------------------------------------------------------------*/ /** @internal @brief Ensure that the vector has the required number of distinct values @param self The vector to check @param ndistinct The (positive) number of distinct values to require @return CPL_ERROR_NONE iff the check is sucessful */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_vector_ensure_distinct(const cpl_vector * self, int ndistinct) { cpl_ensure_code(self, CPL_ERROR_NULL_INPUT); if (ndistinct > 1) { cpl_vector * tmp; const double * dx; double xmin; const int n = cpl_vector_get_size(self); int i, j; if (ndistinct > n) { return cpl_error_set_message_macro(cpl_func, CPL_ERROR_DATA_NOT_FOUND, __FILE__, __LINE__, "A %d-element vector cannot have" " at least %d distinct values", n, ndistinct); } assert( n > 1 ); tmp = cpl_vector_duplicate(self); (void)cpl_vector_sort(tmp, 1); dx = cpl_vector_get_data(tmp); xmin = dx[0]; i = j = 1; do { if (dx[i] > xmin) { xmin = dx[i]; j++; } } while (j < ndistinct && ++i < n); cpl_vector_delete(tmp); if (j < ndistinct) { return cpl_error_set_message_macro(cpl_func, CPL_ERROR_DATA_NOT_FOUND, __FILE__, __LINE__, "%d-element vector must have " "at least %d (not %d) distinct " "values", n, ndistinct, j); } } else { cpl_ensure_code(ndistinct > 0, CPL_ERROR_ILLEGAL_INPUT); } return CPL_ERROR_NONE; } /*----------------------------------------------------------------------------*/ /** @brief Copy a rectangular memory area from one place to another @param self Preallocated location to copy to @param src Location to copy from @param size Size of each element [bytes] @param nx Number of elements in x direction in source @param ny Number of elements in y direction in source @param llx Lower left x position (starting with 1) @param lly Lower left y position (starting with 1) @param urx Upper right x position @param ury Upper right y position @return CPL_ERROR_NONE if OK, otherwise the relevant #_cpl_error_code_. @note self must have place for (urx-llx+1) * (ury-lly+1) elements of size size @see mempcy() Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT if an input pointer is NULL - CPL_ERROR_ILLEGAL_INPUT if the window coordinates are not valid */ /*----------------------------------------------------------------------------*/ cpl_error_code cpl_tools_copy_window(void *self, const void *src, size_t size, int nx, int ny, int llx, int lly, int urx, int ury) { /* FIXME: Need to do pointer arithmetic */ const char * csrc = (const char*)src; cpl_ensure_code(self, CPL_ERROR_NULL_INPUT); cpl_ensure_code(src, CPL_ERROR_NULL_INPUT); cpl_ensure_code(size != 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(llx <= urx, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(lly > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(lly <= ury, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(ury <= ny, CPL_ERROR_ILLEGAL_INPUT); /* assert(sizeof(char) == 1); */ if (llx == 1 && urx == nx) { /* Number of bytes in one row of self */ const size_t rowsize = size * (size_t)nx; (void)memcpy(self, csrc + rowsize * (size_t)(lly-1), rowsize * (size_t)(ury-lly+1)); } else { /* FIXME: Need to do pointer arithmetic */ char * cself = (char*)self; /* Number of bytes in one row of self */ const size_t rowsize = size * (size_t)(urx-llx+1); /* Byte just after last byte to write to self */ const char * cstop = cself + rowsize * (size_t)(ury-lly+1); cpl_ensure_code(llx > 0, CPL_ERROR_ILLEGAL_INPUT); cpl_ensure_code(urx <= nx, CPL_ERROR_ILLEGAL_INPUT); /* Point to first byte to read */ csrc += size * (size_t)((llx - 1) + (lly-1) * nx); for (; cself < cstop; cself += rowsize, csrc += size * (size_t)nx) { (void)memcpy(cself, csrc, rowsize); } } return CPL_ERROR_NONE; } /**@}*/ #define CPL_PIX_SORT(a,b) { if ((a)>(b)) CPL_DOUBLE_SWAP((a),(b)); } /*----------------------------------------------------------------------------*/ /** @internal @brief Optimized median computation for 9 values arrays @param p array to sort @return the requested median @see cpl_tools_get_median_double() Formula from: XILINX XCELL magazine, vol. 23 by John L. Smith The input array is modified. In case of error, the #_cpl_error_code_ code is set, and the returned value is undefined. Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT */ /*----------------------------------------------------------------------------*/ static double cpl_tools_get_median_9double(double * p) { cpl_ensure(p, CPL_ERROR_NULL_INPUT, 0.00); CPL_PIX_SORT(p[1], p[2]); CPL_PIX_SORT(p[4], p[5]); CPL_PIX_SORT(p[7],p[8]); CPL_PIX_SORT(p[0], p[1]); CPL_PIX_SORT(p[3], p[4]); CPL_PIX_SORT(p[6],p[7]); CPL_PIX_SORT(p[1], p[2]); CPL_PIX_SORT(p[4], p[5]); CPL_PIX_SORT(p[7],p[8]); CPL_PIX_SORT(p[0], p[3]); CPL_PIX_SORT(p[5], p[8]); CPL_PIX_SORT(p[4],p[7]); CPL_PIX_SORT(p[3], p[6]); CPL_PIX_SORT(p[1], p[4]); CPL_PIX_SORT(p[2],p[5]); CPL_PIX_SORT(p[4], p[7]); CPL_PIX_SORT(p[4], p[2]); CPL_PIX_SORT(p[6],p[4]); CPL_PIX_SORT(p[4], p[2]); return(p[4]); } #undef CPL_PIX_SORT #define CPL_PIX_SORT(a,b) { if ((a)>(b)) CPL_FLOAT_SWAP((a),(b)); } /*----------------------------------------------------------------------------*/ /** @internal @brief Optimized median computation for 9 values arrays @param p array to sort @return the requested median @see cpl_tools_get_median_9double() Possible #_cpl_error_code_ set in this function: - CPL_ERROR_NULL_INPUT */ /*----------------------------------------------------------------------------*/ static float cpl_tools_get_median_9float(float * p) { cpl_ensure(p, CPL_ERROR_NULL_INPUT, 0.00); CPL_PIX_SORT(p[1], p[2]); CPL_PIX_SORT(p[4], p[5]); CPL_PIX_SORT(p[7],p[8]); CPL_PIX_SORT(p[0], p[1]); CPL_PIX_SORT(p[3], p[4]); CPL_PIX_SORT(p[6],p[7]); CPL_PIX_SORT(p[1], p[2]); CPL_PIX_SORT(p[4], p[5]); CPL_PIX_SORT(p[7],p[8]); CPL_PIX_SORT(p[0], p[3]); CPL_PIX_SORT(p[5], p[8]); CPL_PIX_SORT(p[4],p[7]); CPL_PIX_SORT(p[3], p[6]); CPL_PIX_SORT(p[1], p[4]); CPL_PIX_SORT(p[2],p[5]); CPL_PIX_SORT(p[4], p[7]); CPL_PIX_SORT(p[4], p[2]); CPL_PIX_SORT(p[6],p[4]); CPL_PIX_SORT(p[4], p[2]); return(p[4]); } #undef CPL_PIX_SORT #undef CPL_DOUBLE_SWAP #undef CPL_FLOAT_SWAP #undef CPL_INT_SWAP